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1. Introduction. It has been known for a long time that the lubrication of solid 
surfaces by fluids is primarily an effect of their viscosity. Today, however, there is a 
widespread belief that there are other, as yet unknown, properties of liquids entering 
into the picture, particularly at high pressures. An investigation of such effects is severely 
hampered by the fact that there are few experimental results that can be checked against 
theory. This is due partly to the difficulties of making precise experimental measure- 
ments, and partly to the insufficient development of the conventional theory. 

This paper applies conventional theory to a few situations of experimental interest. 
In Sec. 2, a comparison is made between the Reynolds equation and the more accurate 
Stokes equations. Section 3 studies the cylindrical bearing without assuming film 
lubrication, while Sec. 4 studies the spherical bearing as a sample of a case with side 
leakage. 

2. The equations of Reynolds and Stokes. A good starting point for the study of 
lubrication are the Stokes equations of hydrodynamics,’ 


Vp = nV’V, (1) 
with their immediate consequence 


Vp . (3) 


The relations (1) to (3) are already a special case of the general hydrodynamic 
theory. They assume: (a) that the inertia of the liquid can be neglected; (b) that the 
liquid is incompressible; and (c) that the viscosity 7 is a constant independent of applied 


stresses. 

Assumption (a) is probably justified if the equations are used primarily for the study 
of stresses, for one can easily verify that the pressures obtained from Bernoulli’s theorem 
are of an order of magnitude smaller than the pressures significant in lubrication theory. 
This still leaves the possibility that the flow pattern in regions of low stress may come 
out wrong because of this neglect. Assumption (b), on the other hand, is very serious 
indeed. The Eqs. (1) and (2) always lead to negative pressures that may equal thousands 
of atmospheres. We know with a high degree of certainty that such pressures cannot 
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1Lamb, Hydrodynamics, Cambridge University Press, 1932, p. 595. 
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exist in an actual lubricant.” The most plausible way by which the fluid presumably 
escapes these pressures is cavitation. But a coherent hydrodynamical method which 
takes account of cavitation is not available at this time. We are limited, therefore, to 
qualitative guesses concerning such a development. Assumption (c) is also incorrect 
since it has been shown that viscosity is dependent on pressure and shearing rate. It 
seems reasonable at this time, however, to expect no qualitative change from a more 
refined treatment of viscosity. 

Almost all the calculations on lubrication found in the literature proceed not from 
the equations of Stokes given above, but from a two-dimensional equation of Reynolds.* 
Reynolds’ treatment is known to be more approximate than the treatment of Stokes 
and to be applicable mainly to thin lubricating films; but beyond this, their mutual 
relation is not well understood. 

In order to clear up this point we shall consider the following specialized situation. 
Suppose we have a thin layer of liquid which is bounded on one side by an ideal sta- 
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Fic. 1. A sample case for examining the Reynolds equation. 


tionary plane, which we take as the zry-plane (Fig. 1), and on the other side by an arbi- 
trary surface 


z = h(z, y). (4) 


Let this surface be moving with a velocity whose components are U, V, W which are 
functions of xz and y. 

*Swift, Proc. Inst. Civ. Eng. 233, 267 (1932); Scott, Shoemaker, Tanner & Wendel, J. Chem. Phys. 
16, 495 (1948); Briggs, Johnson, & Mason, J. Acoust. Soc. Amer. 19, 664 (1947); Fisher, J. Appl. Phys. 
19, 1062 (1948). 

3. Reynolds, Phil. Trans. Roy. Soc. 177, 157 (1886). 
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We now introduce power series developments. For the pressure p we set 


1 ~ ; 
7 Pit, Ys 2) = > zx, yz’, (5a) 


v=0 


and for the components u, v, w within the liquid 


u(x, y,2) = DI E(z, ye’, (5b) 
v(x, y,2) = Do a(a, ye’, (5c) 
w(x, y,2) = Dy s(x, ye”. (5d) 
Equations (1), (2) and (3) then become 
Or, dt, , dH, 
dx Ox” T oy” + + 20 + Dé» (6a) 
on, Tm 4 THY + H+ 1) (6b) 
dy —s- Ox” oy” F 4 wees 
ge 3 : 
(+ Varer = SS + The 4 2 + Does (6c) 
OX oy 
\y dg, On, _ 
+ Iie + ar + ay = 0, (6d) 
Or, On, 
b 2 + a i (v + 2)(v + l)z,+2 = @. (6e) 
Ox Oy 


The boundary conditions are 


to = mo = So = 0, (7a) 
> E(x, y)h’(x, y) = Ula, y), (7b) 
Yo ale, Mie, ) = Ve w), (70) 
D fe, Wh, v) = We, v). (74) 


In-reducing the system (6) we observe that the quantities ¢, can be eliminated and 
that Eq. (6c) is redundant, except when v = 0. This leaves the reduced system 
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On, _ Ft, , It, 9 
Ox = Ox" + ay” + (v + 2)(v + 1)é, +2) (8a) 
dn, _ 9m , In, 9)( . 
oy ~~ ax? + ay” + (v + 2)" + 1) n,+2 ’ (8b) 
oF TL 4 DO + Devs = 0 8 
92? ay? y+ 2)v )tr+2 = O, (8c) 
— _ % _ 9m 
i dy’ (9) 
with the reduced boundary conditions 
dX E(x, yh'(x, y) = Ula, y), (10a) 
> n(x, yh'(x, y) = V(x, y), (10b) 
1 fat, 2a} le. ) = —Wle 
p> : { ae oS ay h(x, y) = —W(a, y). (10e) 


We can look at these equations in the following way: the Eqs. (8) determine all the 
£, » and = of index larger than one in terms of 7 , 7, , £1 , m . These four quantities 
are determined by the remaining four equations. An actual procedure of this type is 
generally difficult because (10) involves not only these four variables, but all higher 
index variables as well. 

This is the point where the thin film assumption enters. The Eqs. (8) give us in- 
equalities of the type 


‘ h? 
+2 2 
es Well Mee, ? h” Kax,h” ~ 2,2” 


provided 


; «3. 
Here h is the thickness and 1 the length of the liquid film. This means that in each 
parity, higher powers in the series (5) are small compared with lower powers, and we 
can get a first approximation to a solution by retaining the first term of each parity in 
(10). This leads to the simplified boundary conditions 


U, (11a) 


fh + &h° 


V, (11b) 


I 


mh + neh’ 


1 fat, , an\i. 1 {8 only = 
1 {9% + on, + 3 we +) h* = —W. (11c) 
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They have to be combined with the first of the Eqs. (8a) and (8b): 


OT — 
or = 2b» (11d) 
Oro _ 


We thus have the five equations (11) in the five unknowns &, , & , m , "2, To - If 
we eliminate the first four unknowns from the set, we get 


LJ (, do) , 9 ()39m0)\ _ ow _ poh _ yah (au ov) 
6 ‘2 () are) 4 oy (i oy )} - a9 U Ox ’ wt" Ox + 3y F 


The left-hand side is already Reynolds’ equation. For the right-hand side we can dis- 
tinguish essentially two cases. 

(a) We may have a stationary curved surface rotating against our fixed plane 
(Fig. 2). 
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Fig. 2. Situation for the normal Reynolds equation. 
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Fia. 3. Situation for the reversed Reynolds equation. 
(b) We may have a bump in the moving surface moving parallel to the plane (Fig. 3). 
In the first case we have that 


oh oh 
T= U—+V— 
" U Ox ) Oy’ 


and the last two terms are smaller by a factor h/a, where a is the radius of rotation. 
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In the second case W and the last two terms vanish. Thus we observe that we get the 
Reynolds equation either way, but in the second case the sign, and hence the location 
of the pressure maximum and minimum, are reversed, as compared with the usual case. 
In an actual arrangement h is likely to vary for either reason. We can then make the 


definitions: 


h, = a geometric h defined for an ideally smooth surface; (12a) 
h, = the remainder of h due to roughness. (12b) 


We can distinguish h, and h, by the fact that the velocity of the moving wall is by 
definition parallel to its geometric surface, so that h, is stationary, h, is variable. 


From these considerations, we deduce the final equation: 


ah (n° ors) +2 (n° ons) «oe 4 7k. (13) 
6 \dx Ox oy Oy J Ox OY 
This is Reynolds’ equation if h, = 0. In addition to deriving this equation we have 


obtained a complete description of the flow and pressure pattern through (11) and (9) 
and the recurrence relations (8) for higher terms. This refinement can be applied, for 
instance, to the Sommerfeld solution for the cylindrical bearing or to the results of 
this paper on the spherical bearing.* The result is that the approximation proceeds in 
powers of h/a, where h is the clearance and a the radius, and that there are always 
pairs of contributions, such as £, and £, , which are of equal order. 

It is clear that the restrictive assumptions made in deriving (13) are unsatisfactory. 
If one investigates curved oil films, one gets the obvious generalizations of (11) and (13): 


V= U+ U.N, (14a) 
] 
, Vmo = 202, (14b) 
7 
: V-(h° Vn.) = V-Vi(h, — h,), (15) 
) 


where all vector quantities are understood to be tangential to the liquid film. The 
attempt to remove the restriction illustrated in Fig. 1, however, was not successful. 
But it stands to reason that (14) and (15) should be derivable even if both surfaces 
have irregular shape and proper velocities. The work of later sections and general 
plausibility arguments indicate that the proper generalization of (15) is 
: V-(h°Vn) = (Vi + V.)-{¥h, _ ae vi}. (16) 
We have thus been able to show, at least for a special situation, that the Reynolds 
equation arises from the Stokes equations in first approximation if all quantities en- 
tering the Stokes equations are expanded in powers of the film thickness. In addition, 
we have obtained simple formulas which give such unknowns as the transverse pressure 
gradient from the Reynolds solution.* The standard assumption that the latter is zero 
is, therefore, superfluous. 


4See Eqs. (40) and (59). 
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3. The cylindrical bearing. It has been pointed out by Duffing and by Reissner® 
that the equations of Stokes can be solved rigorously for the infinite cylindrical bearing. 
The derivation to follow will proceed to the end, giving explicit expressions for the 
resultant forces and torques. Several limiting cases will be presented, one of which is 
the classical Sommerfeld solution.° 

The infinite cylindrical bearing has flow in two dimensions only (no side leakage). 
The Stokes equations (1) and (2) reduce, for such a case, to 


ap _ (= oy) - 
ax "Naz? * ay)? (178) 
ap _ (2, #2) 7 
dy Nox + dy" /’ (17b) 

Ou Ov 

oo are ae 

Ox «Oy s (18) 


for the interspace between two circles of radius a, and a, . The circles are placed ec- 
centrically by an amount e (Fig. 4). At each boundary the normal velocity must vanish 
and the tangential velocity must equal a specified constant amount v, and v, , respectively. 

The equations are solved most conveniently by first solving (18) identically through 
the introduction of a stream function V: 


ov 
“= ay? (19a) 


ov 
det te (19b) 


Then if we substitute (19) into (17), we find that we have to solve the equation 
V’vV'v =0 (20) 


and that p will be obtained from VY’W by the method of conjugate functions: 
; a , 
VvV+i »P = f(x + iy). (21) 


The conditions under which Eq. (20) must be solved are that: (a) Y must be a con- 
stant on each of the two bounding circles; (b) the normal gradient of ¥ on each specified 
circle must equal a specified constant, v, or v2 ; (c) the conjugate function of V*¥ must 
be single-valued over the entire domain. 

Conditions (a) and (b) are somewhat weaker than the standard boundary conditions 
for (20).’ Condition (c) which replaced the missing boundary conditions might be 
expected to be awkward. Actually, it gives no trouble in the following applications. 
5G. Duffing, Z. angew Math. Mech. 4, 296 (1924); H. Reissner, Z. angew Math. Mech. 15, 81 (1935). 


6A. Sommerfeld, Z. Math. Phys. 50, 97 (1904). 
7See Frank-v. Mises., Differentialgleichungen der Physik, vol. I, pp. 845-862. 
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i?) 


The solution of (20) may be constructed with the help of a general theorem’ which 
states that any solution of (20) can be written in either of the following two forms: 








Y= 9%, + %, (22a) 
or 
Vv = (** + y')®, + ©, (22b) 
Fig. 4. Basic data for the cylindrical bearing. 
where 


V’e, = V’*?, = 0 (22c) 
and zx and y are Cartesian coordinates. This theorem associates (20) very closely with 
the Laplace equation 


VV’ = 0. 


The electrostatic solution for the geometry of Fig. 4 is well-known. In order to express 
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it simply we pick a coordinate system as shown in Fig. 5. The y-axis is laid through the 
centers C, and C, of the two circles to an external origin O. The distances OC, and OC, 
will be denoted, respectively, by d, and d, . The origin is located in such a way that 


dj —ai = d; —a3 = 8". (23) 
Equation (23) and the definition 
d, _— d, = €é (24) 
give us d, , d, , s in terms of the original quantities a, , a, , e of Fig. 4. We find that 
] bx a a ae (25a) 
5 = De \% 1 9% a 
d, = 2 (a; — ai) + Es (25b) 
a ae : . 
’ = i (a2 — a, — e)(a, — a, + e)(a, + a, + e)(a, + a, — e). (26) 


The points A and B on the y-axis, whose distance from 0 is s, serve as origins for 
the logarithmic potential which is adapted to the present geometry: 


> 
= Cle - (27) 


The curves on which ® is constant are circles, two of which are the bounding circles 
discussed above. If we introduce the variable parameter a for the radius of such a circle 
and 6 for the distance of its center from the origin, we then get the generalization of (23), 


& —a’ = 8’, (28) 


and 
ESS . (27b) 


lw 
&= 5 Cle 


This solution ® and its derivatives can be used in conjunction with the theorem (22) 
to produce a large number of solutions of (20). The successful combination is found to 


have the form 





x +(s+y) y(s + y) y(s — y) 
vY=A]! . aon eee. B= +; == . 3 
°8 x + (s — y) +83 +(s+y) TOR +@e-9 
(29a) 
2 , 2 
+ Dy + E(x? + y? + 8°) + Fy log ,teig, : 
or, more briefly, in a (4, y) coordinate system, 
—— oi +8 sty 8—y 
y= A logs + Base tae at 
(29b) 


6+s 
6-—s 





+ E-2iy + Fy log 
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Here A, B, C, D, E, F are constants to be determined. The structure of the solution is 
clear from (29b); we have picked all solutions of (20) which are linear in y on either 
limit circle. One solution of this form was not admitted, however, namely 





(x* + y° + 8°) log espe = 2y6 log ; + : 


a a 








a 
x 











Fic. 5. Geometric parameters for the eccentric bearing. 


One verifies easily that a term of this type gives rise to a pressure term which is not 
single-valued. It must, therefore, be rejected according to condition (c). 
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The six constants A, B, C, D, E, F are determined by the boundary conditions (a) 
and (b). Since W is linear in y on either limit circle, we get from (a) two conditions, 
annulling the linear part on either circle. Now d¥/dn is found to be linear in 1/y; this 
gives rise to four more conditions, two of them setting the coefficient of 1/y equal to 
zero and two prescribing the value of the constant. Thus we end up with the following 
six equations: 


d, +s 


























1 1 1 1 
Siac" ~8E-0° FO Te (30a) 
ee Se dy +8 py _ 
24482 cee TOT te eS (30b) 
ld-sp late . 
2A +oa pee tad, mad = 0, (30c) 
ld, —s ld, +8 _ 
mito e eae” le sanied 
ld, —s _ild+s Pad 2 pe 
7 es B ca =. y 2a,E + 2sF = —ay,, (30e) 
1d, = = — ld, +. 8 = 2 a Vil 
oa 33 4 4 2a,E + 2sF = —a,.v,. (30f) 
These equations yield: 
—— ] _— 2 
A= 9 (d,d, — 8) 
' ia 2(d, — d;) (a, + a2) (31a) 
(a; + a;)[(az + a3) log {(d, + 8)(dz — 8)/(d, — 8)(d2 + 8)} — 4se] 
4 a;a3(a; Vv; — G2 v2) } 
s(a; + a3)(d, — d,))’ 
B = (d, + s)(d, + s){same curly bracket}, (31b) 
C = (d, — s)(d, — s){same curly bracket}, (31c) 


_ d, log {(dz + 8)/(d2 — 8)} — dp log {(d, + 8)/(, — 8}, 
D= (a + a’) log {(d, + 8)(d. — ad — s)(d, + 3)} may oo (a,v,; + G22) 





= 2s[(as — a;)/(a> + ai))(aws + 202) (314) 
(a; + a2) log {(d, + s)(d, — 8)/(d, — 8)(d2 + 8)} — 4se 





a;a3(a; ‘Vv; — a2 V2) 


(a; + aye ' 
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a (1/2) ) log | {(d, + s)(d, — 8), fh — s)(d, + 8)}(awi + ar) (316) 
~ (aj + @; ) log {(d, + s)(d, — \/(d, — s)\(dz + 8)} — 4se ’ 


e(aw, + av.) (31f) 











(a; + a) log { iA (d, + s)(d. — s)/(d; — s)(d2 + 8 )} — Ase 


From the flow pattern, which is completely specified by (29) and (31), we can de- 
termine the pressure by using (21). We find that 





prety , (8 — y) _ py __48t (32) 


Ba coe _¢2 | 
nt y(6 + 8) y'(6 — s)° y(o — 8°) 





Calculation of the forces requires, of course, the entire stress tensor. In our case, the 
tensor has three components which are obtained as follows: 














sa satis. 
To = ~?P ” 9a ay’ 
my 
~~ - ae 
Tw idl Ox dy’ 
_ (2% 2%) 
Fon ™ % Ox" ay /]* 
The result is: 
1 1A ree _ - 5 _ Ss) B a[2y’ + (s — 36)y — 288] | 
—T.. = —4, - Se SOS ——— 
n y(e — s*) y(6 + s)° 
(33a) 
— P 286 7 sr 95 aah 2 ,2 
4 a czy — | (s + Soy +285] _ gp 82(2 by - : - 35 +s 2 
y (6 — s)° y(s — 8°) 
1 sa(2iy — & — 8°) a[2y’ + (3s — 5)y + 2s") 
=, = 44 SP Ss 5 pee a 
gee yea y(5 + 8)° 
(33b) 
wT) ini 8 2s e7(25 -_ 2 3” 
_¢ a[2y” "— (88 + 6 y + 2s] 4 4P5 x(2 Dey - i — #) . 
yb — s)° y(e — 8°)? 
1 s[ —25y° + ( (35° r+s Jy — 26s") 
~., = 4A i 3 ll 
n y"(o — 8s) 
ez Pp + 2(26 — s) sy" + (36s — 2s" = 5)y + s(6 — 8) 
y 76+ s)° 
(33¢) 
, —2y* + 2(26 + sy’ — (8és + 28° + Oo )y + 8°(6 + 8) 
y°(6 — s)° 
4 4F§ —2éy’ - + 45°'y — 6(3° + 8°)] + 8°)] | 


y(é — s°) 
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The main interest centers usually on those resultants of (33) which act on the pin 
or bearing, when considered as rigid bodies. The equations for the components of force 
acting upon a circle of parameter 6 are: 


L $ (x. = + Wry v=) dl; ’ 
a a 


I 


+F, 


L § (x. Sim u—1) dl, - 
Qa 


+F, 
a 
We find for these expressions: 
+F, = 8rnLF = W, (34a) 
where F is given by (31f) and L is the length of the cylindrical bearing. Similarly, 
F, = 0. (34b) 


This equation is due to the anti-symmetry of the pressure pattern about the line of 
centers (cf. Fig. 10). This same symmetry precludes the transport of z momentum 
into any of the circles, and thus makes the force independent of the parameter 6. The 
sign depends, of course, on the side from which the force is acting. 

Substituting (31f) into (34) we get explicitly that 


a 8arnLe(aw; + ayv2) F (35) 
(a? + a3) log {(hh + 8)(dz — 9d, — 9d, + 8) — 400 





W = 


The geometrical parameters entering into this equation are explained in Fig. 4, Fig. 5 
and Eqs. (23) to (26). 

The total torque depends on the choice of the fulcrum. Taking the fulcrum at the 
center of the circle of parameter 6, we find that 


2 — — 2 
+1 = Lat {(re — x,) SY 42,7 a. } dl, , 





or that 
+T = 8xnL(A + FB), (36a) 


where A and F are constants given by (31a) and (31f). This reduces for the inner surface 
to 


i = 8arnL(A + Fd,) (36b) 
and for the outer surface to 


T, = —8enl(A + Fa,). (36c) 


If torques are referred to the same fulcrum throughout, such as the origin, then a constant 
torque 


+T = 8rnLA 


results. This constancy can again be checked from symmetry considerations. 
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A convenient way of representing the resultant torque is by calculating the co- 
efficient of friction f. The generally used definition is 
j= td (37) 
: a,;W ’ 
where 7’; is the torque referred to the center of either circle and a; is its radius. It should 
be emphasized here that we get two coefficients of friction, one for the inner cylinder 
and one for the outer cylinder. If we assume the outer surface stationary, we find for 
the inner surface that 
al(dids — s*){(a? + a2) log {(d, + 9)(ds — 8)/(d, — 8)(dy + 8)} — 400 














a 2a,80"(a? + a) 
(38a) 
__ ae 
a; + a; 
and for the outer surface that 
a,(d,d. — s°)[(aj + a>) log {(d, + s)(d, — s)/(d, — s)(d. + s)} — 4se] 
fh. = - , a ps 3 
2se"(a; + a2) 
(38b) 
~ a tay 


The f generally calculated in the literature is f, , although many experimental arrange- 
ments actually measure the torque on the outer surface. 

A good deal of insight into the nature of the solution can be gained by considering 
limiting cases. We are starting off with the most important one. 

(a) The Sommerfeld-Reynolds limit: e = d, — d, — 0. 
We expand our results in powers of e and retain only the first term. Let us define the 


nominal clearance ¢ as 


Q2ao~-&=C (39) 
and write further that 

a, ~a,™~ a. 
From (25) 

d~d~= 
and from (26) 

_ 4r2_ 21/2 

e"- (c "Sei 
Then we find for the pressure from (32), 
_ Gye(aw, + av2)(2c — e cos ) sin y (40) 





’ 


alee (2c? + e’)(c — e cos y)” 


where y is the azimuth measured along the oil film in a counter-clockwise direction, 
starting from the point of minimum clearance. 
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Equation (40) is obtainable directly from the Reynolds equation (11) and has been 


so obtained by Sommerfeld.° 
The same limiting process can be carried out on the resultant force and torque. 


For the total force we get 


- 124nLea’(v, + v2) 
= Qe +eVe —e)” (41) 





W 


with the usual ambiguity in the sign. For the torques we find that 








_ 4rnLa’{v(c? — e) — v,(c* + 2e*)} 
T, = (Qe? > ee ane ey ? (42a) 
2 2 2 . 2 2 
T, = 4nnLa {v,(c — e) — v(c + 2e°)} , (42b) 


(2c? + e’)(c pea oe 


d= a2 











Fia. 6. Limit of zero clearance. 


The coefficients of friction reduce to 


fi 3 (43a) 





lc —e 
Pa =3 ‘ (43b) 
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The formulas for W, 7, and f, can be found in standard textbooks. The flow pattern 
reduces in this limit to a superposition of viscous drag by the inner cylinder plus a 
Poiseuille type flow from the high to the low pressure region. This latter flow is opposed 
to the motion of the pin on the side of wide clearance. 

(b) Limit of zero clearance: s — 0. 

This limit, shown in Fig. 6, has some interest. From the experimental side, we observe 





Fic. 7. Streamlines for a journal bearing of diameter ratio 5:9. Case of no eccentricity. 


that a loaded bearing usually comes very close to a situation where there is contact 
along a line or in a point. Mathematically, the simplification obtained is considerable, 
mainly because the basic solution (27) is now rational instead of logarithmic. In addi- 
tion, the distinction between d; and a; , 6 and a, or c and e disappears. As a consequence, 
the entire system of equations (29) and (31) is now replaced by the short expression 


V=- ?.. {a,v,(@ — dz) + AV2( a) }(@ ne a;)(@ = a2) (44) 
ac 
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and (32) is replaced by 
_ _ 2naia; x (2 2,2 2 ) M1 ws} 
y= e ay {o + V2) y - + ds + a, a, de » (45) 


This case is very favorable for a study of flow conditions. The infinities in the stresses, 
on the other hand, give rise to an infinite force and infinite inner torque. The coefficient 
of friction stays finite, however, and is 


f, = 3S a =. (46) 





Fig. 8. Streamlines for a journal bearing of diameter ratio 5:9. 50% eccentricity. R = points of flow 
reversal on stationary surface. 


Another finite result for this limiting process is 
T, = 0. (47) 


Neither (46) nor (47) has much practical significance, since each'is obtained by a 
cancellation of infinities. 
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(c) Limit of concentric cylinders: d, ~ d, ~~; d, — d, — 0; a, , dy finite. 
This complicated passage to the limit gives well-known results. We find that 
1 asd, — ids 2, a102(0 "rs + avs) 


‘a -e a = an log r, (48) 





Fic. 9. Streamlines for a journal bearing of diameter ratio 5:9. Full eccentricity (zero clearance). 


where r is the conventional radius. The pressure and load are, of course, zero. There 
remains a shearing stress in the r, ¥ direction which is 


¢ ye -1 1 ‘ 
2naj;a,(a; v; + as Vo) 1 
‘y= ee aime 


=. 49) 
a; — a; r* (49) 
27rL? times this expression is the exerted torque. 

(d) Limit of a plane and a cylinder: dz , dg >; 8, d; , a finite. 
This interesting limiting case does not alter the nature of the solution obtained, but 
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greatly simplifies all coefficients. We shall, therefore, keep Eq. (29) and replace (31). 
Writing 











d, = d, a, = a, vy = 2B, vy, = V, 
we get 

dV 1 dav 
4¢-Dgitau-a ts" (50a) 

__ _ 2@d+s)V (d + s)av 
asia log (d + s)/(d — s) s ‘ (50b) 

a 2(d — s)V ; (d — s)av 
* log (d + s)/(d — s) 8 : (50c) 
D= -YV, (50d) 
E =0, (50e) 

' V 

F = (50f) 


log (d + s)/(d — s) 

The solution decomposes into two parts. The V-part describes a situation where a 
cylinder falls with velocity V near a vertical wall; the v-part applies when a stationary 
cylinder rotates with speed v. 

If we substitute the expressions (50) into (34a) and (36b) we observe that (a) a 
cylinder falling in the neighborhood of a vertical wall experiences no torque from it, 
and (b) a cylinder rotating in the neighborhood of a wall experiences no force from it. 
In case (a) there is a retarding force equal to 





ae 8arnLV 
~ log (d + s)/(d — s) (51) 


In case (b) there is a frictional torque equal to 


W 


T = 4anL dav ; (52) 


s 


As an illustration of the formulas obtained, numerical results are shown for a bearing 
having a diameter ratio 5:9. Figures 7, 8 and 9 show the flow pattern for three different 
eccentricities: Fig. 7 shows the concentric case (Eq. (48)), Fig. 8 a case of intermediate 
eccentricity (Eqs. (29) and (31)), and Fig. 9 the case of zero clearance (Eq. (44)). The 
salient feature of these pictures is the return flow going in a direction contrary to the 
moving pin. In Fig. 9 the entire flow returns that way, whereas in Fig. 8 only a fraction 
does. This flow is due, of course, to the pressure difference which forces liquid through 
the wide section as well as through the narrow clearance. Even in the Reynolds limit 
of film lubrication this type of flow persists. 

Figure 10 shows the pressure distribution for the flow condition of Fig. 8. We can 
verify on this graph that the pressure generally varies little in direct line from one surface 
to the other. The Reynolds treatment neglects such variation altogether; this assumption 
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is thus fairly well justified, even for a bearing of fairly wide clearance. The anti-sym- 
metry of the pressure pattern is more doubtful since it is directly connected with the 


assumption of incompressibility. 





| | 
+4 = es 
4 42 0 -e 4 
Fic. 10. Lines of constant pressure for the bearing of Fig. 8. Pressures are to an arbitrary scale and 
contain an additive constant. 


Finally, in Figs. 11 and 12, are shown load vs. eccentricity curves and friction vs. 
load curves for this bearing (Eqs. (35) and (38)). In order to compare them with the 
traditional Sommerfeld formulas, an intermediate “radius” of 7 is ascribed to the bearing, 
and Sommerfeld’s equations are applied. The differences are not very marked. The 
graphs do emphasize, however, that there is a difference in the coefficient of friction 
according to whether it is defined on the inner or outer circle. 

4. The spherical bearing. When we pass from rotating cylinders to rotating spheres, 
we pass essentially from a bearing with line contact to one with point contact. As a 
consequence, we get “side leakage” of the liquid about the contact point. The pre- 
sentation to follow gives a partial solution of the mathematical problem. The difficulties 
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encountered are not surprising since the solution for cylinders is based so much on the 
potential (27) for the same geometry. 

In discussing such a bearing there are essentially two possibilities: either the axis 
of rotation is parallel to the line of centers, or it is perpendicular. Bearings of the former 
type have been constructed and examined recently by Shaw and Strang.” While the 
bearing apparently works very well, there is at present no complete understanding of 
its operation. It is easy to show from Stokes’ equations (Eqs. (1) and (2)) that if the 
liquid motion is entirely azimuthal, the pressure in the crescent-shaped cavity is a con- 
stant. Actually, there is a slight upward flow in the bearing which is force fed. It may be 


2 oe a See 


Ww 
e) 














| nar 





Fig. 12. Coefficient of friction versus load for a journal bearing of diameter ratio 5:9. The left graph 
defines friction as usual on the inner circle; The right graph defines it on the outer circle. 
- Present calculation — —-—-— Sommerfeld formula 


that here we have simply a tightly sealed-off pocket of liquid which is under the hydro- 
static pressure introduced by force feeding the oil. This would make the bearing hydro- 
static rather than hydrodynamic, except perhaps for a hydrodynamic effect which pre- 
vents contact at the lateral rim. 

The normal mode of loading a bearing is lateral. Let the z-axis be the line of centers, 
and let the z-axis be parallel to the axis of rotation (Fig. 13). The angular velocity will 
be denoted by w. The remainder of the nomenclature will be the same as in Fig. 5. 

If we introduce a cylindrical coordinate system, 


p=2r+y, (53) 


y - 
tan ¢ = " ; (54) 





8M. C. Shaw and C. D. Strang, The hydrosphere, a new hydrodynamic bearing, J. Appl. Mech. 15, 
137 (1948). P ‘ 
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then the fundamental equations (1) and (2) read, 





_— Vv», — 4, -37°. 
=? - V0, — “56 + 35" 





a 


Fic. 13. Coordinate system for the spherical bearing. 


on the stationary sphere: p? + 2° — 2d,z + s’ = 0, 


v, = 0, 
v, = 0, 
v, = 0; 


on the moving sphere: p’ + 2° — 2d,z + s° = 0, 
v, = —w(z — d,) sing, 
vl. = w(z — d,) cos g, 


Vv, = wpsing. 
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The form of the equations and boundary conditions is such that the dependence on 
gy can be eliminated from the equations. We find that 


v, = f(p, z) sin ¢, (55a) 
Vv, = g(p, 2) cos ¢, (55b) 
* v, = h(p, z) sin ¢, (55c) 
p = q(p, 2) sin ¢. (55d) 


AXIS OF ROTATION 





EQUATORIAL 





PLANE 








Fic. 14. Polar coordinates for the Sommerfeld type calculation of the spherical bearing. S = point of 
smallest clearance, L = point of largest clearance. 


At this point difficulties are encountered. It is relatively easy to eliminate g and gq, 
but the reduction to one variable does not seem possible. This is the place where the 
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introduction of the stream function solved the corresponding cylindrical problem. An 
analogous trick does not seem possible in this case. However, the symmetry of the 
flow pattern is sufficient to permit the prediction that the resultant force is at right 
angles to the axis of rotation and the line of centers. 

Since the general problem appears insoluble, we shall solve it in the Reynolds limit. 
This also has not been done up to this time. We place a spherical polar coordinate 
system into the center of the two spheres, pointing with the polar axis toward the point 
of least clearance (Fig. 14). 

For such a coordinate system the components of the boundary velocity are 


V> = —wa sin ¢, (56a) 
V, = —wa cos 8 COs ¢. (56b) 
Substituting these expressions into (15), we get 
[ #i,. op . ot, se a 
—— = ¢h> sin 3! = = = —G6na’w sin g = - 
sin 8 a8 {h im? sa) * aber \” ae CMomes, (57) 


In the coordinate arrangement of Fig. 14 h is only a function of 8; in good approxi- 


mation, 


h =c — ecosd, (58) 


where c is defined in (39) and e in (24) or Fig. 5. Equation (58) shows, in agreement 
with (55), that p has the form 


p(d, ¢) = qd) sin ¢, 


where q(#) satisfies an ordinary differential equation. One can verify by direct inspection 
that this equation has a simple particular integral which is 
3na*we (c — e cos 8/2) sind 


qv) = e+ e*/4 t~e cos 8)" 


This integral is the only one which is finite at d = 0 and vd = r. It is, therefore, the correct 
solution for the closed spherical bearing. The pressure distribution in such a bearing 
is thus 
6nea*w(2c — e cos 8) sind sing 
_ oned =. (59) 


Poe) = “ad + Ole — e cus 8)" 





There is a remarkable similarity between (59) and (40), particularly if we specialize 
(59) to give us p along the equatorial plane (see Fig. 14). In this plane ¢g = +7/2, and 
the azimuth y essentially equals 8. More precisely, we must substitute 

cos ? = cos y, 
sin 3 sin g = —sin y. 
Using these substitutions, we find that the dependence of p on the azimuth y is exactly 
the same in the two cases. There is a minor difference in magnitude which shows up 
best by forming ratios: 
p(59) _ 2c* +e (60) 
p(40) 4¢e +e 
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In view of the inequality 


Gs 24< 4 
this means that 
59) ; 
a <2. (61) 


ae ~ 
2 = p(40) ~ 5 


It is clear from this analysis that side leakage does not play a very essential role in 
the theory of lubrication. All the main features are the same in the two cases: the 
pressure maximum near the point of least clearance, the region of negative pressure on 
the exit side, the anti-symmetry of the pattern. The factor sin ¢ gives the continuous 
transition from the entrance to the exit side, and a factor of the order 1/2 indicates 
the effect of side leakage on the pressure. 

The symmetry of the problem specifies, as pointed out earlier, that the resultant 
force is at right angles to the line of centers and the axis of rotation. In magnitude, the 
contribution of 7) , as given by (59), outweighs all others; this may be verified by a 
study of Sec. 2 or by direct computation.” The resultant W is obtained from the formula 


p2e 


W = | sin 3 dé | dy p(?, v) sin 3 sin ¢ 


“0 0 


which yields 





6rna‘w f 2 9 ct+e 
W = 3 (CC + €) log —— — 2e 62 
e’(4c’ + e*) | +) bce ‘ (62) 
where a, c and e are given in (39) and Fig. 4. Now a) makes no contribution to the re- 
sultant torque. We must pass, therefore, to the next order, which is due to the shearing 
stress (terms for wu which are linear and quadratic in z; see Eq. (5b)). This frictional 
torque is parallel to the axis of rotation and is obtained on either surface from the 
formula 
: Ro i ‘aad ' av, Ov. 
+T = na® [ sin 3 di | dy {sin y + cos? cos ¢ ive} 
“ r r 


“0 0 


The derivatives in question are obtained from (14), (56) and (59). The result is, 
for the inner surface: 


4 2 2 : , 
7, « Ste G +7) eg tt - zee (63a) 
e€ Ac’ + ¢ c—e 
for the outer surface: 
4 2 9 
an wna w 2(2c° e°) f » < c+e \ : 
a a ee a y ; a: ae HO 6: I 
Ts 3 47 +e G + e°) log pene 2cey (63b) 


As above, (Eq. (42)), the torques are not exactly equal and opposite because the ful- 
crum is in the center of each sphere. If we reduce to the same fulcrum, with the help 
of (62), the torques check. 





*Bull. Nat. Res. Council 84, 1932, p. 236. This publication calculates the effect of shear for the 
Sommerfeld case, but ascribes the wrong direction to it. The more rigorous treatment is also very inap- 


propriate as compared with our Sec. 3. 
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The coefficients of friction become 


_ 2c’ +e’) 

fi . 3ae ’ (64a) 
Qc? — e* 

<< “* (64b) 


There are some practical difficulties in making the flow pattern of this bearing 
visible. The general nature of the flow, as given by (14), is the same as for the Sommer- 
feld bearing; there is a viscous drag flow which varies linearly across the film. This 





Fig. 15b. Lines of current flow for the spherical bearing. 
A = axis of rotation, S = point of smallest clearance. 


drag has a tendency to accumulate material on one side; the accumulation is removed 
by a Poiseuille type flow which varies quadratically across the film. This flow is more 
complex in the present case than in the Sommerfeld case, for it can go at any angle 
to the movement of the bearing. The best over-all picture of the flow is gained by showing 
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the total stream flow Q, that is, the velocity integrated across the film. It results from 
the formula 


Q0,¢ = | VO, ¢,2) a 


or from (14a), 


Q = 3U,h’ + 3U,h’. (65) 





Fic. 15b. Lines of current flow for the spherical bearing. 


A = axis of rotation, L = point of largest clearance. 


Q is a divergence-free vector in two dimensions and is, therefore, derivable from a stream 
. . . . . * 3y/ 
function VY which measures the flow between two points on the spherical film in cm*/sec. 
And ¥ is obtained from (60) by either one of the two equations 
av a 


= —aQ, or 


— = aQ; sind. 
ad F) Qos 


8 14 
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If we apply Eqs. (14b), (14a), (65) and these last relations successively to (59), we 
find for VW: 

_ __ awe 

~ 4¢ +e 
Two photographs of a plot of this flow pattern on a sphere are shown in Fig. 15. The 
case picked is the one having the most extreme side leakage when e = c. The side leakage 
around the point of smallest clearance (marked S) is clearly visible. It is also seen that 
the flow does not take place about the points A which mark the axis of rotation but 
about a point which is displaced to the side of large clearance. The vortex points, ob- 
tained from (66) by setting 





5 (2c — ce cosd — e’) sin 8 cos ¢. (66) 


max 
Vd "te 
(@, ¢) ale, 
are 
le , 
COS? max = — dc F SIN ~ mx = 0. 


This indicates a maximum displacement of the vortex by 30°. The total flow between 
vortices comes out to be 
rs (c? Ss e’ /4)°” 
’ (68) 


Fe _ Y via =a  ¢ + e’/4 








The flow is thus always smaller than it would be without eccentricity, but there is no 
change in order of magnitude, in contrast with the cylindrical case. 

It should be emphasized that Fig. 15 shows only the integrated flow across the liquid 
film. The actual velocity varies widely from layer to layer. As an example, we see in 
Fig. 16 the flow in the layer adjoining the stationary sphere. The flow lines obey the 
equation 


Pe sin 3 COS ¢ — (69a) 


13c7/(c?+e?) 


3 2 2 3 
—2c’ — ce cosd + 4ceé — € cosd | 





or for the special case e ~ c, 3d K 2/2, 


F COs 





The point of minimum clearance is marked M in the figure. There is little resemblance 
between Fig. 15 and Fig. 16. The essential feature of the latter are the points R which 
separate two regions of behavior. Along the line RMR the direction of fluid motion 
coincides with the direction of bearing motion, but outside of R along the same straight 
line the direction is reversed. We can observe this phenomenon for the cylindrical 
bearing in Fig. 8, where the corresponding points are also marked R. It is obvious from 
Fig. 8 that an integrated flow pattern cannot give these essential details. 

5. Conclusion. This paper has given formulas of greater rigor and generality than were 
available heretofore for the cylindrical bearing. As an example, we may mention the 
new formula for the load carried (Eq. (35)) and the new formulas for the coefficient of 
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friction (Eqs. (38)). The spherical bearing is brought up to the status which the cy- 
lindrical bearing had up to this time by the Eqs. (62) and (64). It is unfortunate, how- 
ever, that these results are of limited usefulness because of the negative pressure difficulty. 

The approach of the research worker to this basic problem has been rather sterile. 
The negative pressure region has been referred to as “‘inoperative’’’’, and the general 


———— ll \ 
Paci = 


/ 7 - 





} 

a Se en 

Fic. 16. Flow pattern in the spherical bearing. Flow direction at the stationary surface. Region limited 

to point of minimum clearance. D = direction of moving bearing surface, M = point of minimum clear- 
ance, R = points of flow reversal. 








practice has been followed of simply ignoring this region in all integrations. This pro- 
cedure yields essentially a frictional torque divided by 2 and a load divided by 2 sin a, 
where a is the angle between load and line of centers. The procedure is intuitively 
reasonable, but is open to the objection that it is not mechanically sound. Any cut-off 
line leaves unbalanced stresses acting across it; and if the line of centers is chosen for 
this purpose, there remains an unbalanced shearing stress which is responsible for the 
increase in the load obtained. This seems to indicate that we must first understand 
better the manner in which the liquid cavitates in the “inoperative” region. Then the 
hydrodynamic equations can be reexamined and corrected accordingly. New values for 
force and torque will thus result which will be mechanically consistent. 

In the absence of such a theory we must adopt a cautious attitude in applying 


0H. W. Swift, Proc. Inst. Civ. Eng. 233, 267 (1932). 
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theoretical formulas. It seems plausible to assume that the order of magnitude of the 
pressure maximum and of its distance from the line of centers would be reasonably safe 
against modifications because the differential equation is the same on the high pressure 
side of the bearing. Limiting ourselves to the Sommerfeld limit (although the general 
case can be handled too by observing that the F term in (27) is generally small compared 
with the other two and that B and C are very simply related), we can calculate the 
maximum of (40) and its location. We find that 





COS Wimax = re 3 (70a) 


or 





1 _ [Qc — ce — 2" 
tan 5 vex = E: Folero] _— 

Formula (70) applies to all loads; but if we are mainly interested in high loads, we 
can take the combined cases, limit (a) and limit (b) only: 


c~we, h=c-eXe. 


We then get for the distance D of the pressure maximum from the line of centers, 


a v - (2 ay 
D 2a tan 5 = 3c a. 


This is written in a more useful form from (39) as 


iw (2 wea) (71) 


It is very satisfactory that only the difference in curvatures enters into (71). This makes 
D dependent on local conditions only, as postulated. If we substitute (70) into (40) 
and make the same approximation as above, we get 


= 27 1/2 nv 
Pus = (24) (1/a, — 1/a,)'"h*” (72) 


This formula also depends on local conditions only, and thus we may expect (71) and 
(72) to hold under much wider conditions than the derivation given would indicate. 

A simple alternative derivation which does not use the film lubrication condition 
proceeds, for instance, from the plane and cylinder case, limit (d). From (32) and (50) 
we find that the pressure along the plane has the following form: 





1 x 
rn p = —8avrd @+e). 
Assuming as before that d ~ a > s, we get from (26), 
2 2h 
8 la ~-Ua" 


Maximizing p and substituting s, (71) and (72) follow exactly. 
Calculations proceeding from more general assumptions have been prepared by the 
author. Thus formulas (71) and (72) have so wide a sweep that they can very probably 
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be applied to convex types of contacts where the curvatures would add, of course. More 
attention must be paid in this case, however, to other stresses which become more 
nearly comparable, as the results of limit (d) show. 

The difficulties which the theory encounters in the field of lubrication indicate 
strongly that both theory and experiment should concentrate part of their attention 
on the negative pressure problem and its possible connection with cavitation. A clear 
mechanical picture of hydrodynamic lubrication, backed up by quantitative measure- 
ments, is probably essential, before modifications of this picture can be discussed with 
hope of success. 

Acknowledgment. In conclusion I wish to thank Mr. Murray L. Deutsch and Mr. 
Thomas A. Lindland for their valuable help in carrying out the computations. 


STUDY OF THE BOUNDARY LAYER AT SUPERSONIC SPEEDS IN TURBULENT 
FLOW: CASE OF FLOW ALONG A FLAT PLATE* 


BY 
CARLO FERRARI 
Cornell Aeronautical Laboratory 


INTRODUCTION 


The flow of a gaseous fluid past a rigid, smooth half-plane, parallel to the asymptotic 
velocity of the stream, is studied. The edge of the half-plane (leading edge) is assumed 
to be orthogonal to the stream, and, in addition, it is assumed that the flow is every- 
where ‘‘turbulent’”’. The motion within the boundary layer is investigated in an attempt 
to clarify many of the assumed and conjectured facts relating to the turbulent layer 
within a supersonic stream. For instance, it is found for supersonic flow, within the 
range of Mach numbers examined, that a drag law holds which is formally identical 
with the one which is true for incompressible fluids, and, moreover, that the thickness 
of the boundary layer also obeys a law which is formally analogous to the one describing 
the growth of the layer in the incompressible fluid case. In addition, such relationships 
as those concerned with the velocity and temperature distributions within the layer, 
and the transfer of heat from the stream to the constraining plate at various Reynolds 
and Mach numbers are deduced from clearly set forth premises. Analogies and extensions 
are presented for full understanding of the phenomena described. 

In presenting the results of this study of the motion within the turbulent layer it 
is found necessary to obtain the solution as an asymptotic one, i.e., one which is rigorously 
valid at a sufficient distance from the leading edge. The formulas obtained contain a 
parameter, the determination of which must be left to experiment. A plausible value of 
it, based upon the results obtained by Fréssel’ from his investigation on sub- and super- 
sonic streams in pipes, will be given. 


GENERAL EQUATIONS 


1. Symbols. Bars and primes will be used to denote, respectively, the mean values 
of the following quantities, and their fluctuations about their mean values due to 
turbulence: 7 = enthalpy; 7) = enthalpy of the external stream at zero speed; 7, = enthalpy 
in contact with the wall; 7, = enthalpy of the external undisturbed stream. 

In the next set of symbols, the subscripts 0, p, a, and the bars and primes give the 
quantities the same meanings as those specified above for the enthalpy: 7 = absolute 
temperature; U = component of the velocity along the z-axis; V = component of the 








*Received Dec. 20, 1948. Italian manuscript submitted to Cornell Aeronautical Laboratory July 
22, 1948. This work was carried out under Navy BuOrd sponsorship. 

IW. Fréssel, Strémung in glatten geraden Rohren mit Ueber- und Unterschallgeschwindigkeit, Forsch. 
auf dem Geb. des Ing.-Wesens 7, 75-85 (1936). 
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velocity along the y-axis; W = (U? + V’)'”; p = density of the fluid; 1 = viscosity 
coefficient; y = kinematic viscosity. 

In addition, the following symbols are used: p = pressure; r = tangential friction 
stress; r, = tangential friction stress at the wall; U* = friction velocity [U* = (r,/p,)'”"); 


V, = limit velocity of the undisturbed stream; u, = U./V; ;f= U/U*; f. = U,/U"*; 
o = U/U, = f/f. ; u* = U*/V, ; M, = Mach number of the undisturbed stream; 
y = T./T, ; 6 = thickness of the boundary layer; t = time; c, = coefficient of local 


friction drag (c, = 2r,/p,U2); C, = coefficient of total friction drag for a plate of length 
x; q = amount of heat transferred through the plate from the fluid, per unit length of 
wall. 

The meaning of other symbols used in this paper will be given as they appear in 
the text. 

2. Equations of motion and of continuity. In the region outside the laminar sublayer 
one may neglect the transport of momentum due to thermic molecular vibration in 
comparison with that due to turbulence. Let the z-axis coincide with the line of inter- 
section of the rigid plate with the plane on which the motion is being studied. Then 


one may write: 


op — a, r) J T?) < TV) 
_ ay vl ) + ax (pU") + ay (pUV), 
(1) 
Pp _ 9) wi 2 (yy) + & (7? 
- ay (P! )+ ar (pUV) + ay (pV°). 


Letp =Dt+p,p=pt+p,U =U+ UU’, and V = V + V’. Then at a sufficiently 
large distance from the leading edge of the plate it is true, in all likelihood, that p = p(y), 
U = U(y) and V = 0. Therefore, if one neglects, as is usually done, pV’ in comparison 
with pUV, it is found that 


pUV = pU'V’ + Up'V" + p'U'V’ = +r = const. = 7,, 


p = const. = Dp. . 
It certainly seems allowable to admit that p U'V" may be neglected in comparison 
with the other terms; then one can substitute the simplification 
pU’ Vy’ + Up’ Vy’ = T, (2’) 


for the first of Eqs. (2). 
From the equation of continuity, 


Op fe) 
ot - Ox 


(oU) + 5, (eV) = 0, (3) 
the following relation may be established: 
pV’ = cp,U*, (3) 
where c is a parameter, related to the Mach number and concerning which we know 
solely that 
lim c = 0 as M,-— 0. (3’’) 
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3. Energy equation. If in addition to neglecting p’ U’ V’ one neglects, as in Sec. 2, the 
viscosity and the molecular conductivity in the region outside of the laminar sublayer, 
one easily may find that 


aT @ (2) -1(% , Op éP) 
Cat Z(e ae a eee Y oe” (4) 


where c, denotes the specific heat at constant volume. From this equation we obtain 


1 (i+w) =2. 


Pat ~ Ot 
But 
a oe ee ie )] af (i 1 )| 
pS (i+ gm )- 2 | pu{i+ 3 tay Leni te® 
ol (i42 ”) 
+ at | o + 2 W ’ 
and, therefore, at a sufficiently large distance from the leading edge, one has 


pV + 3W*) = ~ pV + 4W’)’ + + 4U")p’V’ = const. = ap,U**. = (4’) 


4. Explicit form of the terms defining the turbulent transport of momentum and 
energy. In order to describe the terms which define the turbulent transport of momentum 
and energy as functions of the geometric and kinetic characteristics of the average 
motion, we may go back to the well-known considerations which permit the establish- 
ment of the following expressions for these terms:* 





75 tan dU Y772\1/2 i 1 x | = fj 1 7) yr2y\1/2 
UV? = 1a, Vli+an =lglitgl (Vv. (5) 


where | and 1, are lengths which are likely to be related to the scale of turbulence. It 
appears plausible to put 1 = 1, . In fact if 

i + 3U* = const. (6) 
(the flow is everywhere isoenergetic), the correlation between V’ and (i + W*/2)’ may 


be assumed equal to zero, and therefore V’(i + W?/2)’ = 0. On the other hand the 
following relations hold: 








Va + 3W’) = Vi +4VW? = VV" +-27(20U' +: U" + C&V”) 
ae >77 Tr wr?’ — di J7/2\1/2 Thy’ 2\'/2 dU , 
=wV7i +UV'l adr” yy’? + UV") > (6’) 
_ eed, 4 1 =). 
= (: dy 2! dy 0, 


2G. I. Taylor, The transport of vorticity and heat through fluids in turbulent motion, Proc. Roy. Soc. 
London (A) 135, 685-705 (1932). 

































36 CARLO FERRARI [Vol. VIII, No. 1 





where the meaning of the symbols is obvious. Now, in order that Eq. (6’) be a necessary 
consequence of Eq. (6), whatever the particular variation of 7 and U in the region may 
be, it suffices that 1, = 1 everywhere, and therefore also that 1, = L. 


VELOCITY AND TEMPERATURE WITHIN THE BOUNDARY LAYER 


5. Law of variation of temperature within the boundary layer. Taking into account 
Eq. (3’) and Eq. (5), Eqs. (2’) and (4’) become 


, dU pay 1/2 Tk _ = *2 
a (V + cp, U ppU ? 
(7) 
a =. a +z rj + (i re v*) pp,U* = ap,U* 
p dy 1 5 & a 1 5) cp,U* = ap, ; 
One obtains from them that 
i+ 3U* =i, + BU, (8) 


where 7, is the enthalpy at the wall, while B is related to the parameters a and c through 
the equation 





B+ cy — aU* = 0. (9) 
la at the edge of the boundary layer, where i =i, U = U,, one must have i, + 
U2/2 » + BU, . It follows immediately that 
p= kt Te—§ a3 _ i (,_1) 
U, U, U, Y 


(10) 


II 


ee 5 +e DM 
where k is the adiabatic exponent and y = T,/T, . The law of variation of temperature 
in the inside of the boundary layer is therefore obtained as 

















T _1_U*_(k—1)Ms_ r=1/1 1 |Z 2(k — 1)M?2_ 
Ty a+e-pMet 27 k— DMG. 2+ — DM: 
| - (11) 
1 YU {r= - (k —1)M?2_ U 
ae v yao tee | 


This relation is identical with that defining the law of variation of 7 within the boundary 
layer in the case of laminar flow when the Prandtl number is equal to unity’, Pr = 1. 
6. Law of variation of velocity within the boundary layer. 
6.1. Integration of the equation of motion. In order to carry on this discussion further 
the quantity 1(V’’)'’”, denoting the so-called coefficient of turbulent transport, must 


8Th. v. Karm4n and H. S. Tsien, Boundary layer in compressible fluids, J. Aero. Sci. 5, 228 (1938). 
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be given an explicit form. Now the same formulation will be assumed for this coefficient 
that is commonly accepted in the case of incompressible fluid flow (one might repeat 
here the reasons which justify it).* 

Let, therefore, 

UV")? = Uty®(f, f’, +++) = Utyeo + af + af’ + +++), (12) 
where f = U/U* while f’ = df/dn and ap, a, a, °°: area priori unknown functions 
of » = y/4, if 6 is the thickness of the boundary layer. 

Accepting the same degree of approximation as in the case of incompressible fluids, 
let 


UV™)'? = U*yay , (12’) 
where a is a constant. Hence 
¥772\1/2 dU _ Tr*2.. ¢7 
lV’) rs = U*'nf'a. (13) 


Now, since it is true that 














bbb __[,_1B Be Ze) 
> %T i, —U*/2+BU 2U; i tp U. \i, 
(14) 
_ _ l a ——$ $$$ 
1 — uyo + (y — Io’ 
the first of Eqs. (7) becomes 
{p,U*?/[1 — uzyo”? + (y — Vel}aonf’ + cp,U**f = p,U*’ (15) 
or also 
f’ *2 £2 J ss *2 £2 Fi 
aonf’ + cff 1 — yu" f + f GY- Di =1— wf +F (y — 1). 
By separation of variables one obtains 
ay df dn 
sea . =—. 16 
(lL — yu® f+ —-DS/f) —cf) 1 (16) 
Let 
mo ( _ )"] 
f= yu, E =i + fy — i) Us ’ 
(17) 
_— <t1| ( dy y"] 
fo — Qy'/" Ua 1 + 1 + (y le 1)’ Ua ’ 


where one employs either the upper or the lower sign before the radical according to 
whether y > 1 ory < 1. Then we have 





4G. Moretti, La teoria del trasporto turbolento, L’Aerotecnica 27, 44-50 (1947), and Considerazioni 
sullo studio della turbolenza nei tubi e nei canali, L’ Aerotecnica 27, 404-410 (1947); or C. Ferrari, Lo stato 
attuale della teoria della turbolenza, Lectures on Phys. and Math. in Univ. of Torino, 1936, p. 223. 
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I ae a 
(fi — 17 u® fer” 


u* f + fo(1 — cf) 


[1 = yur’f? + Gy — DF/F0 — cf) © 


sis nite en: _— 
fy — 7 u* f + yutt + fe at ef’ 
where 
yu ut 
j= he SS 
in? fa(y'u* == Ca) i (hi + fay’ u* + cf.) ’ 

—¢’ (18) 

Se hr 

(y'"u* + cf2)(y'"u* — cf) 


One should observe that within the boundary layer f/f, < 1, while f,/f.y'’u* > 1, 
—f./fay'°u* < 0, for any value of y whatever; thus the following inequalities certainly 
hold: —f./fay'’*u* < f/fe < f:i/fay'u*. One has, besides, that f/f, ¥ 1/cf, ; and this 
relation, which may now be accepted as a restriction which the constant c must satisfy, 
will be justified later in Sec. 10.2. 

Thus it appears that the function of f appearing in the first member of Eq. (16) has 
no singularities within the layer. By integrating Eq. (16) one obtains the equation 
a| - = ; log (f, — y'7u*f) + of 5 log (y'*u* f + fr) — 3 log (1 — of | 

yu yu c 
(19) 
= log 7 + log A, 


where log A denotes a suitable constant of integration. 

6.2. Determination of the constant of integration. Equation (19) describes the law of 
variation of the mean velocity in the inside of the turbulent boundary layer, and, there- 
fore, for all values of f included within the interval, f, < f < f. , where f, denotes the 
value of f at the boundary of the so-called laminar sublayer.’ For values of f close to 
f. , for which one has y'u*f Kf, and y'"u*f <f, , it is permissible to put 


log (f; — 7’ ‘u* f = log fi — yur , log (y’ ‘u*f + fr) = log fz +7 w=, 


log (1 — cf) = —cf. 
Taking into account Eqs. (18), one may put Eq. (19) in the form 


ail? \ ' f \ 51/2 : ; 
=e | -(x u* + f.)( log fi — 7" "a" = ) — (e u* — /)(Iog fo + 7" """ r) 
Cc Cc , Jiu Cc ; ; Je 


(19’) 
“] mp l/2 1/2 
— (fi + fe cf | = (f, + fy(? u* + f.)(? u* — f, (log n + log A). 
- 


At the boundary of the sublayer, the thickness of which will be denoted by y, , one has, 


therefore, 


’Th. v. Karman, Turbulence and skin friction, J. Aero. Sci. 1, 1-20 (1934). 
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M1Y whAth, rae h 
u lati 


f? 
c hi fo ° . 


7a meal ef; + ak 


vel Lipteatnan: 
it fe 14 + fu)(¥" “ - f,)(log % + log A). 


And as f,f, = 1, there is obtained also: 


afi + jo] u*? + i u*(f, — fo) — |p 
res ° or. + (v'22 — 5) to 5 
“* Ye 
(fi + fe) Cae a f(y” a: — f,)(log r} + log 4) 


= (f, + die u** + y¥ v2 ue > MSs —fi)- 1|(t08 % + log 4), 


whence 


* * 
f.= (og ve + log A) a | -(r <4 fn) log fi + (1 ae 
Qo \ 0 e 4 Cc 


c 


en ee ee oer) 
f.) log 5 | (y7u* /c f.)(y*u* /e */c — f,) Qo log ** 5 + log A (20) 


mn Vom, Map See a ee y.U* 
~ fi + fs) ( yut/e — fi log fi + x ut/e + fe log f,) = ~~" 


where it is assumed, as is done by K4rman’, that within the sublayer the velocity is a 
linear function of y. Thus the following equation may be derived: 





, a 0 l ; log fo 
log A = —log : + a = es Ti a0 ( ay + th _), (21) 


When Eq. (21) is inserted into Eq. (19), we find the relation 


r log (1 — 7'u* f) 





Ee ae 1 a } 1/2, * ) 

(fi 4 fo)(y'*u* Cc fr) log (1 ite fo (22) 
sania inital 1 ee snails ane * 

+ Gu" /e + fo u*/e — fa log (1 of) | = log y* + B 
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where ‘ 


yU* y,U* ‘ , 
y+ = and B= a & with y, = (22) 
1 P 
Vy Vy Pp 


It is assumed here, as in the study of incompressible fluid flow, that a) and y,U*/», are 
universal constants, which are independent, in particular, of the value of y. It should 
be observed that for y = 1 one has f, = f, = 1. Thus Eq. (19’) leads, for the value of 
log A given by Eq. (21), to the following fact: 


f= (log y* + B). (19’’) 


Ao 


This relation is identical with that holding for incompressible fluid flows. This must 
be true not only because the equations established in the case of incompressibility must 
be deducible as a particular case from those which are true for gases, but also because 
for y = 1, in proximity to the walls, the effect of compressibility is negligible. 

Thus it seems allowable to assume in every case, for a as well as for B, the values 


obtained from experiment for p = const. 


DRAG LAW AND LAW OF VARIATION OF THE THICKNESS OF 
THE BOUNDARY LAYER 


7. Dependence of the local drag coefficient upon the local Reynolds number. 
Equation (22) gives, as already said, the law of variation of the velocity outside of and 
on the boundary of the laminar sublayer, i.e., for y > y,. . By following a procedure 
identical with that used by Prandtl” it is, however, possible to extend the applicability 
of this law right up to the wall and thus evaluate the coefficient of local resistance, 
c, = 27,/p,U, . To this purpose, a transformation of coordinates is performed such that 
for y* = 0, f = 0. One need only substitute log (1 + by*) for log y* + B. When this is 
done, it may easily be verified that for suitable values of 6 and for y* sufficiently large 
the new expression is practically coincident with the old one, and that the new expression 


vanishes for y* = 0. Equation (22) thus becomes 
ao - . ] = - 1/2, x r) 
Cc | (fi + fe)(y'*u*/e — fi) log (: Ts if 
SS es a == log (1 + y*u* r) (22’) 
(fi + Jody “a /e + J2) fo 
re log (1 p | log (1 + by*) 
Ey Sar ee Te ae —<¢ = lo y™). 
(y' 2u*/eo + fo (yu e— f;) Ig CJ) Og \ IY 


*L., Prandtl, ““The mechanics of viscous fluids,” Aerodynamic Theory, ed. by W. F. Durand, vol. 3 


Julius Springer, Berlin, 1934, p. 145. 
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Equation (22’) may be used in order to express the coefficient c, in terms of the Reynolds 
number, Rk, = U,2/v, , by means of 


-4f 00 7 
ae pU(U, — U) dy. (23) 
The method explained hereafter is similar to that used by Prandtl.° Remembering that 


t, = p,(U./f2), dy = (v,/U*) dy* and putting 6* = U*s/v, , we find from Eq. (23) that 


- oes - d > ~~ = m= &€ Fs dy* ¥ 
Pepi = Poe ge | 101-0) dy = p07 | ® fa — 2) Fede. (28 


: *v0 


Now from Eq. (22’) one obtains 


] + a ae 2u®/e+fs) 
ee Ce es 
] + Y ( ais yu* f/f, 
(24) 
x (; —_ yu* f/f ; 
where 
ft=fhth. (25) 
Let 
U* u 
ai Saad | 26 
c= 8 V, 8 yf.’ (26) 


where 8 is a numerical parameter, a priori unknown. It should be noted that the re- 
lationship (26) represents merely a heuristic hypothesis, which is the most simple repre- 
sentation that may be formulated in connection with condition (3’); but nothing justifies 
a priori the assumption made in this paper of regarding 8 as constant. Equation (24) 
can be written in the form 


1/2 ; fa/uef* 1/246f,) 

l+y¥ eal fs)" sila tie 
pa a. on: 

l1— y“ou,/f; 


_ Bu. @ Bf a/ualy*/*+BSs) (y*/2*—Bhs) 
< (15,0ue_) | 


1+ by* = ( 
(24’) 


172 


— you/f; 
By differentiating Eq. (22’) one obtains 


1/2 1/2 
I y“u* 1 » el he 





b dy* “ ut ghd, = 
1+ by* c Face ec— ffi — yy rutf © ft(y'Put/e + fe) y'utf + fr 





(y'"u*/e + fa)(y'*u*/e — fi) 1 — ef * 


ol ay df 
(fi — Vout fy ut f + fol — ef) 





_ sit f, do : 
~ [1 — yuio® + (y — Ieol(1 — Bue) ’ 
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so that taking Eq. (24’) into account, it is found that 


dy* Ay f = l 7 : a 
do b** {1 — yuo + (y — Lel(1l — Buc) 
(27) 
1+ y" “TU, Ta ewe 1 — Buco a 
Xx wn oS 7a 77 , 
l—y¥you,/f; L— yy" ot./ Ss 
in which 
a, = a/ftu(y'? + Bho), A = aoB/uly'” + Bf)(y'” — Bf,). (28) 
By inserting into Eq. (23’) the expression for p/p, given by Eq. (14), one obtains 
_ a d f2 [ a(1 — oa) 1 
vf. bdzt** Jo [1 — ucyo” a ty = l)o]* (1 — $c) 
(29) 
: l +y¥ ou, A ge 1 — Bu,o ners 
x = 173 ip | .. een do 5 
\l — you,/ J; l= ota Ii 
By differentiation it is found that 
U, Gi. [ o(1 — a) (; +- 7" “OU, a) re 
<= ar - Fz - = = : 
et P b ay" [1 — usyo” + (y — Lo}(1 — Bue) \1l — y'ou,/ fi, 
f 1 — Buse — ; 1+ vy ‘ou,/f 
> ( 73 2 -) 2+ f.a, log s Ee 
\l — y“‘ou./f; oy otha 7 
(29’) 
fa | — Bu,o } Ho» df, (f.) 
+ ts 10 —— ae Do( Ja); 
ee 1 — 7" “OUa/ Ji co 
whence 
[ da ¥ a2 » ay e * ‘ 
dR, = = fro0(f.) df, = ofa) dfs, (30) 
V; ) ) 
where 
$3 = fGo (30’) 
Finally, by integrating Eq. (30), it evolves that 
¢ a : 5 Qo ‘ 
R.=] olf.) df. = x(f.). (31) 
b Jo ’ i b . 
Observe now that 
[ ( 2p,\' 
i, = =- = ) = [2y(1 — u;)/c,]' 
‘Dis (r,/p Cr Pa 
(32) 
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and let R, = U,x/v, . Then we find that 
— lt+n 
R, = RB, = RB, — ( : ) (33) 
Va 7" l—u, 


which, in conjunction with the first of Eqs. (32) and Eq. (31), establishes the relationship 
being sought between the Reynolds number of the undisturbed stream and the coefficient 
of local friction. 

8. Relationship between the coefficient of total friction and the Reynolds number. 
The total friction drag for a length x of the flat surface, is given by 


F, = | t,a2 = | p,U** dz = p, | a 7 ; fief.) df. 


J0 0 /0 . 


(34) 
= py Us | ilf.) dha = tps Wf.), 
D Je b 
where 
afa 
1 = fo and y = | ¢:(fa) df, . (34’) 
Consequently, one deduces that 
QF 7 ¢s - 
— Fe - 2F,, os 9 ¥ fa) (35) 
p,l at uy UR, x(Ja) 
or also that 
2F j 2 2 ~, 
C= = 2¢ adi ~ OC. = Si - Hf. (35’) 
pal me 2 Pa x 


9. Law of variation of the thickness of the boundary layer. From Eq. (22’), keeping 
in mind Eq. (26), one obtains for y = 6, 


i 
4] eer = :( Vays Le) 
ao} ft u*(y'? — Bf,) we (1 — f,/ + frtu® y'”% + Bho ss fo 


{a 


/ 


I 


; — pee 
un(¥ : + Bfo){(y'”? B) — I] log a as 5 | 


/ e * 
= log ( ] b bl 4 = log (1 oa o Rs), 
\ y fa 


p 


where we have set R; = 6U, v, . It is thus found that 





b ‘T+ yu, Lay" 1 — Bu, y" : 
= us. Sarr eT. ard , 36 
I+ ‘. Rs (j — "Us f,/ l—y"u./fi (36) 
and as, in addition, one has ; 
= U 1 1 , 
By 3 SR = —— (36’) 
Vo y "(1 — u,) 
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Eqs. (36) and (36’) give the parametric representation of the line constituting the edge 


of the boundary layer. 
10. Explicit formulation of the drag law and of the law of variation of the thickness of 


the boundary layer. 
10.1. Numerical results. In order to be able to give an explicit formulation of the 


drag law, one must, first of all, evaluate the integral defining the function ¢(/,). 

















TABLE I. 
M, = 1.5 Ua = 0.5571 B = 0.1 
P Cr ¥o ¢: dR,/de, 
7.744 0.02300 7.83 60.64 — 5.274X10' 
10.81 0.01181 26.37 285.0 —13.07X10° 
13.87 0.007168 92.37 1281.0 2.055 X 10" 
16.94 0.004809 332.1 5624.0 — 2.446108 
20.00 0.003448 1216.0 24320.0 — 2.42410" 
23.06 0.002593 4514.0 104100.0 2.539 x 10" 
M, = 1.5 la = 0.55708 B = 0.25 
fa Cr #0 ¥1 dR, de, 
16.52 0.005055 380.5 6285.0 - 2.413108 
20.00 0.003448 1887 .0 37740.0 - 3.762 XX 10° 
23.48 0.00250 9501.0 223110.0 — 4.97310" 
M, =1.5 Ue = 0.55708 B= -l 
fo Cr o ¢1 dR,/de, 
20.00 0.003448 140.0 2801.0 — 2.798 x 10° 
M, = 1.5 ug = 0 55708 B =0 
fa Cr $0 ¥1 dR,/dc, 
20.00 0.003448 934.0 18680 .0 —18.620 108 
M, =2 Ua = 0.66667 B=0.1 
Te Cy ¥o ¢1 dR,/de, 
16.33 0.004165 199.2 8154.0 — 5.343108 
20.00 0.002778 2681.0 53620 .0 — 9.672X10° 
23.67 0.001984 14710.0 348100 .0 — 1.457X10" 
M, =2 Ua = 0.66667 B = 0.25 
fe Cr 0 ¥1 dR,/dce, 
15.83 0.004432 844.0 13360 .0 — 7.496 x 108 
20.00 0.002773 6960 .0 139200 .0 —25.11 10° 
24.17 0.001902 59450 .0 1437000 .0 — 6.667 X10" 
M,=2 Ua = 0.66667 B=0 
t. Cr Yo 1 dR, /de, 
20.00 0.002773 1875.0 37500 .0 — 6.765 X 10° 


Since it appears impossible 


to work out a closed-form expression suitable for the numerical 





computation of this integral, the following procedure has been followed, which permits 
the establishment, at least qualitatively, of the drag law and of the law of variation of 
the thickness of the boundary layer. Their quantitative determination will be made 
possible only after experimentation is undertaken to arrive at the value of the parameter 
8 introduced above. In order to obtain a rough estimate of such a value 8, reference is 
here made to the above-mentioned experimental results given by Fréssel which deal, 
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however, with friction along tubes rather than along plates. The simple case of a 
thermally insulated wall, y = 1, is now considered. The function g, then takes the form 


— [ ; all a / (3 4. we) "(! Buse)" . 
~ (1 — w%0)(1 — Bu,o)\l — uo l — uo 





“0 


(29’’) 


> |2 + a,f, log 1 + ue + f.a, log Eas but | de. 





1 — u,o 1 — uo 


The values of this function, for different values of 8 and of f, corresponding to 
M, = 1.5 and M, = 2, have been calculated numerically as shown in the Appendix. 
The results obtained are summarized in Table I, where the values of dR,/dc, calculated 




















TABLE IT. 
Me =1.5 B=0.1 
(c,)-1/2 c,°/2 dR, /de, c,5/?(dR,/dc,) ; 
6.594 4.231 3.989 
9.202 19.81 19.91 
11.81 89.39 94.14 
14.42 392.30 397 .00 
17.03 1692.00 1950.00 
19.64 7262.00 8694.00 
M, = 1.5 B = 0.25 
(c,)~W2 c,5/2 dk, /dc, c,5/2 (dR./dc,); 
14.06 438.4 352.8 
17.03 2627 .0 1950.0 
19.99 15560 .0 10852 .0 
M,=2 8=0.1 
c,)~¥2 c,!? dR, /de, c,5/2(dRz/de,) 
15.51 594.8 804.8 
18.97 3933 .0 6030.0 
22.45 25536 .0 44337 .0 
M,=2 B = 0.25 
(c,)~"/2 c,5!? dR,/de, c,5!2(dRz/dc,) ; 
15.02 980.3 615.4 
18.97 10210.0 6030.0 
22.92 105400 .0 58490 .0 
M, = 1.5 B= -1 
(e,)-12 c,5/? dR, /de, 
17.03 195.4 
B=0 
17.03 1300.0 
Moe =2 6 =0 
(c,)~¥2 c,5/2 dR,/dc, 
18.97 2751.0 


with Eq. (30), are also given. These values have been evaluated by taking into account 
the second of Eqs. (32) and by assuming for a» and for b the values 0.4 and 8.93, re- 
spectively (these values are very close to those determined by Prandtl® in his investiga- 
tion). 
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On the basis of the above-tabulated data, the values of (c,)~'”? and of ¢2”(dR,/dc,) 
have been calculated. They are given in Table II where, in addition, the values of 
c,’"(dR,/de,) for the same values of c, , but for an incompressible fluid, are also given. 
These quantities have been calculated by means of the following equation, 


/ IR \ 9! 2 { (Me 1/2 - . "] 
syol alt a ; 1/2 <C,) C. C; an 
e | : — exp {ao(2/c,) 1(1 — +-—)/-—-—], (37) 
Ic. / b Qo a QX 


\c 
which can be obtained as a limiting case from Eq. (30). It is also possible to deduce it 
from the formulas given by Prandtl.” 

The values of (c,)~'’* and of c?” | dR,/de, | tabulated above, are plotted on semi- 
logarithmic scales in Fig. 1. It may be noticed at once not only that such a diagram 
[log (¢ dR,/dce, |); e,'*), referring to an incompressible fluid, approaches very closely a 
straight line within a very large interval, but also that the same situation holds with 
still higher approximation for the diagrams referring to a compressible fluid, and 
equally well for M, l.bBand M. = 2asfor B = 0.1 and p = 0.25. 

Therefore, it seems allowable in all cases to substitute for Eq. (30) the Eq. (38) 
which follows. This equation holds at least for values of the Mach number not differing 


considerably from those given above: 


' 2 | aR, Set 
co?’ =A log (<! F + B, (38) 
\ de 
Note here that A and B are parameters which eventually become functions of the 
Mach number only. From the diagrams in Fig. 1, we infer, by setting A’ = (log,) e)A = 
0.43429 A, that 
M, = 5 8B = 0.1 A’ = 4,1 B = 3.8 A = 1.781 
= 0.25 = 3.85 = 3.8 = 1.671 
M, =2 8B = 0.1 1’ = 4.25 B = 3.7 A = 1.846 
= 0.25 = 3.91 = 3.27 = 1.698 


For an incompressible fluid, in the interval corresponding to the highest values of 


(c,)'’”, one finds A’ = 4.075; B = 3.625; A = 1.7697. 
10.2. Drag law. From Eq. (38) the drag law may be derived at once. In fact, let 
A log B, = B, then Eq. (38) may be written in the form 
> 5/9 dk, 9Q7 
c,'’? = A log (Bue! r ), (38’) 
de, 
whence 
| dR, we 
Be”? | — = exp (1/Ac,’*) 
ac, 
or 
Bee ee 2A ; ‘79 os : 1A” 
.=- ( (i/ Ac.” = exp (1/Ac,’*)(1 — 2Ac,’~ + 2A’c,) — = 
Ms B,J, °xP ate ae Bo 
also, for R, sufficiently large and hence for c, small enough, 
i By ' 
c,'’? =~ A log (R,c,) + A log oA” (39) 
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Therefore it may be stated that: for supersonic fluid flow, at least for Mach numbers 
not too different from those here considered, a law of resistance holds which is formally 
identical with that valid for incompressible fluids (logarithmic law of Karm4n). It 
follows, consequently, that the relationship between the total and the local friction 
coefficients is the same as that given by Karman for incompressible fluids, i.e., 


C, =c¢,(1 + 2Ac'”). (39’) 


A quantitative determination requires the knowledge of the values of A and of B, , 
i.e., of 8. Now the experimental results obtained by Fréssel concerning fluid flow through 
pipes, while confirming that the logarithmic drag law holds also at supersonic speeds, 
give, moreover, exactly the same values of A and B, for Mach numbers larger than 
unity as for the case of incompressible fluids. Assuming this result to be true also for 
the fluid flow past a flat plate, one may carry out the determination of the parameter 
8 from the data given in Fig. 1, where the results obtained for a given value of c, for 
8 = 0and 8 = —1, are also plotted. It follows that the value of 8 for any Mach number, 
at least for those not very different from those considered above, must lie between 0.1 
and 0.25. A further rough estimate of 6 will now be obtained by imposing the condition 
that the various curves, related to different M, , pass through the same point of the 
*/? | dR./de, |); c;'“*] plane and by supposing it to be permissible to interpolate, 


° 


flog (Cc, 






30 





M,= 15, 2 = 
Mg = 15, 8 = 025 
Ma= 2, £ = 0.25 


ab 





x» 





Se 
-C, = 
. dc, 


Fic. 1. Mach number influence in shifting the semi-logarithmic plots connecting the drag coefficient and 
Reynolds number. 
for the values of 6 within the given interval, between the values obtained for 8 = 0.1 


and 6 = 0.25. 
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Thus, by insisting that, for 7, 1.5 andc,'” 17.0286 (corresponding to f, = 20), 
the curve go through the point which is the ordinate of the curve representing the case 
of an incompressible fluid, one obtains by interpolation between the values given in 
Table II that 8 = 0.14. Analogously, by imposing the condition that for ¢,'/” 18.973 
(corresponding to the same value f, 20), the curve related to M, 2 should pass 
through the point of the curve representing the case of 1, = 0 which corresponds to 
the above ordinate, one obtains by interpolation between the values tabulated in Table 
II that 6 = 0.15. 

Figure 2 presents the graphs describing the law of variation of c, ‘7? in terms of the 
log (c.” | dR./de, |) for the Mach numbers M, = 1.5 and M, = 2 and for the value of 
8 = 0.145. They have been obtained by linear interpolation of the values corresponding, 
for the above Mach numbers, to 8 = 0.1 and 6 = 0.25 respectively. 
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Fic. 2. Law of variation of (c,)~/? in terms of log (c,5/2 | dk,/dc, |) for the Mach numbers M, = 1.5 


and M, = 2 and for the value of 8 = 0.145. 


It may easily be verified that: 

1. The diagrams are really very close together and also close to the graph representing 
the case of an incompressible fluid flow (the points of which are denoted on Fig. 2 
by crosses). 

2. The drag coefficient decreases with increasing Mach number. This decrease may 
scarcely be perceived, however, on the semilogarithmic scales, for values of the 
Mach number varying between 1.5 and 2. 

Therefore, the assumption that the value of 8 depends only in very slight measure 
on the Mach number seems to be plausible. This justifies the procedure followed here 
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and warrants the use of expression (26) as well. It should be noted, however, that this 
property is a consequence of the assumption made as to the validity of the experimental 
data given by Fréssel and their extensibility to the case of a flat plate; obviously, 
therefore, it deserves confirmation.* 

10.3. Law of variation of the thickness of the boundary layer. For the purpose of es- 
tablishing an explicit relationship interconnecting the thickness of the boundary layer, 
the coefficient of local resistance, and the corresponding local Reynolds number, calcula- 
tions have been carried out by means of Eq. (36). For the calculations we have used 
M., 1.5 and the values of f, and of 6 already given in Sec. 10.1, which determine the 
values of RP, and of R, c\” . The results obtained are summarized in Table III. In this 
* in the case of an incompressible 


table, for every c, , the corresponding value of R,c}’ 
fluid is also given; it has been calculated by means of the equation 


91/2 


} (exp fa(2/c,)'”? — 1). (40) 
) 


1/2 
Cp R; - 


Taste III, 
M,=1.5 B=0.1 


Cy R; Rs (c,)? (Rs(c,)"/) 5 
0.02300 59.52 9.027 6.472 
0.01181 352.9 38.35 28.72 
0.007168 1893.0 160.30 125.20 
0.004809 9625.0 667 . 50 551.20 
0.003448 17290 .0 2777 .00 2409 .00 
0.002593 226800 .0 11550.0 10600 .00 

M, = 1.5 8B = 0.25 

C R; R;(c,)¥2 
0.05055 9983 .0 709.8 
0.003448 70220 .0 4124.0 
0.002501 478800 .0 23950 .0 


The values of (c,)~'”” and of c:’*R,; are plotted on the semilogarithmic scale in Fig. 3. 
It may easily be seen again that not only in the case of an incompressible fluid, but also, 
with the same approximation and for any value of 8, in the case of gaseous flow it is 


allowable to take 


c,'’? = D log (R; c)’””) + E, (41) 


r 


where D and E denote parameters which are eventually functions of the Mach number 


only. From Fig. 3 one has 


M, = 1.5: B= 0.1 D’! = logye D = 4.17 E = 2.65 D= 181 
= ().25 = 3.805 = 3.19 = 1.65 
M, = 0: D’ = = 4.04 E = 3.39 D = 1.75 


*In the light of recent experimental results for a flat plate the theoretical treatment presented here 
has been improved by the author. The new analysis confirms that there is very appreciable decrease in 
drag coefficient with increase in Mach number, and it is found necessary to select a smaller value of 8, 
near zero, for this flow condition, in contrast to what is found for pipes. 
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Following a procedure analogous to that indicated by Karman,’ the relation being 
sought may be derived from Eq. (41) and Eq. (39). Subtracting Eq. (39) from Eq. (41), 
and observing that the values of D and of A are sensibly equal, therefore putting D = 


A, one obtains 





4; — B’ 2 
6 = exp {_ 2-2 Vell (42) 
( A 
where B’ = A log (B,/2A). Hence, the law of variation of the thickness of the boundary 


layer is also formally analogous to that holding for incompressible fluids. As far as the 


M,=15, £* 0.1 
Ma = 1.5, B= 0.25 





1 0 10* 5 
1G; Ry 
Fie. 3. Boundary layer thickness in relation to drag coefficient 
quantitative determination is concerned, it is interesting to see what results correspond 
to the above values of A, B, E in both cases, i.e., for 8 = 0.1 and B = 0.25 as well as 


for M, = 1.5 and M, = 0 (incompressible fluid). One has 


7 — B’ . 
M, = 0, B’ = 1.388, exp {- BB} = 0.3226, whence 6 = 0.3226zc}’’; 


M,=1.5, 6=0.1, B’ = 1.539, 


exp & _s =F = 0.5354, whence 6 = 0.53542c}’’; 


0° 10* ed 
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M, = 1.5, B = 0.25, B’ = 1.784, 





ul / 
exp {~ -- - a = 0.431, whence 6 = 0.4312c,”’. 
These results show that when @ lies between 0.1 and 0.25, the thickness of the boundary 
layer grows more rapidly than in the case of incompressible fluid flow. Assuming as a 
probable value of 8 the previously given one, 8 = 0.14, and admitting that it is per- 
missible to interpolate, for 8 lying within the above-mentioned interval, between the 
values of 6/xc,’* that we have found, one deduces the following law of variation of the 
thickness of the boundary layer (M, = 1.5): 


= 0.52c'”, (43) 


Note: The results now attained have been obtained assuming y = 1, but from the quali- 
tative point of view they still hold for y ¥ 1. A quantitative determination would require 
that the calculations be repeated, retaining for 8 the value now deduced. 

10.4. Determination of an asymptotic expression for calculating the function ¢o(f.). 
The result obtained in Sec. 10.1 suggests the determination of an asymptotic expression 
(holding, that is, at sufficiently large values of f,), for the calculation of the function 
¢o(f.). In this way the result is justified analytically, and we get a general expression 
which permits us to investigate the influence of the different parameters (Mach number, 
ratio of absolute temperatures, y) on the drag law and on the law of variation of the 
thickness of the boundary layer. Such an asymptotic expression may be obtained as 
follows. 

The integral describing ¢o(f,) may be written in the form 


. +1 7 ai-—e«)  _ |2 ( + 7’ “oute/ fr 
™_ I ( — y'?ou./ fi)’ + y'ou,/ fe)*"(1 — Buse) 2+ falar log — 7"? ou,/f; 
+ a, log a an e. )| exp {1 as log | = i 7 RE 
l — ou./ fy 
aL Q> log 1 —— — ‘ee -_ | do. (44) 
Let 
o=- 1-6 (45) 


one finds that 


1/2 Ua - 1/2 Ua _ 1/2 Ua MA a 
log (1 9 qs “) = log (1 + “)(1 Y fel + yu, )| 


12 Ue) ae se FC ld 
= log (1 9 ts) Y fol + y'Ue/ fe 37; 7 (lL + y'7u/ fo)” 


to 





= log (1 + 7’ : te sal y'” Zi 1 + ayy he gi(e), (46) 
\ J2 2 2 
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where g;(¢) is a function of ¢, which, as ¢ varies within the interval 0 to 1, remains finite, 
continuous, and positive, and which is equal to unity for e = 0. 
Analogously, one can write 


3/2 Uae — " _ 1/2 Ua 1/2 Ug _ a 
log (1 Y yf, c) log (1 Y ts) + Y f, =e yu Tf go(e) , 
(46’) 


log (1 — Bu.c) = log (1 — Bu.) + = egs(e), 


where g2(¢) and g;(e) are functions which satisfy the same conditions as does g,(e). It 
is then true that 


L+¥+ V2 ou, — Bu,o 
exp ae log Teel J t+ sali om /fi 








(47) 
1+ 7'u, Lh)" ( i= Be)" Sale ta 
= (3 poe 1 ul fi 1 oe ual! hi exp { fable} 
where 
YUa/ fo 1 ol! fr _ — mm 
d(e) = 1 +4 i. wy ag, + -_ yu, i i — Bu, A273 (47’) 
and 
Qs; = a, + Qe. 


Within the range of values of 8 given in Sec. 10.2 and for e varying between 0 and 1, 
¢(e) is always positive. 
Furthermore one has 











1+ you, / fo 1 — Bu,o -) 
—* f(a og | you,/f, + ET — you h 
. 1+ y'*u,/f Ll — Pu. ) 
= 2+ fla go ———_- = gs | om Se (ele: 
| flay log — yu if, = Qe log = Ua) [2 f.b(ee; 


=—2 
1/2 Ua ee — - oe 9 1/2 Ua ae am ( \ 48 
(1 ve .) “T+ 77 u/fr OP {2 fo + yu] fa O° ae 
—2 
1/2 Us ee ee 1 “ 1/2 Ua <- 
(1 a fi -) .. toe uy, /RY exp { + fi 1- uel fi ” (a}, 


(1 — Bu,c)* = = exp {- Te casto. 





Thus Eq. (44) becomes 
= ft (3 +7 vateta)" 'G 1 pe)" 
~ (1 — u)?(1 — Buz) \l — yu /S, L— y’"u/ fi 


x [ d — 6 exp [ef + d)]I2 + SF — fucd] de 


0 


(44’) 
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where one has put 








a ee ee ee 
oe er hl — 7 adh Hf) + 7_ pu, 2) — 2 fol + 7" Ua/ fe #9) 
(49) 
ous t+ wl fr __1 — Bu, 
F = a, log Daas Y Ute/ fi a, log 1 Yue fi : 


In order to evaluate the integral in the right-hand member of (44’), the interval 0 to 1 
may be divided into two intervals, 0 to 7 and 7 to 1, where 


n= fe". (50) 


Let H be the maximum value attained by the expression «(1 — ¢)[2 + f,.F — f.ed] 
within the interval 7 to 1. It is found that 


[ e1 — )(2+ f.F — f.eb] exp [—e(f.6 + o)] de < H [ exp [—e(f.@ + ¢1)] de 


(51) 


= — 4 — fexp[-(@ft + af.)] — exp [-f. + 4)]} = OCexp [-of'"), 


fb + ¢, 
which can be neglected for f, sufficiently large. In the intervalO < ¢ < nlet@ = do + Yoe 
and ¢; = ¢1.0 + Wie. We have then that 


exp [—e(f.d + ¢:)] = exp [—(fado + o1,0)€] exp [—(favo + vie]. 


Now, if we denote by y,,,o and by y,,,, the maximum absolute values that the y may 
reach within the above interval, it is found that 


exp [—(fabo + dr.ode] exp [—Vm.ofa® — Varfa”) 
< exp [—e«(f.6 + ¢:)] < exp [—(fabo + d1,0)€] exp [Wnofa*® + Ynirfe’”*] 


or 


exp [—e(f.d + ¢)] = [1 + O(f2'”)] exp [—€(fabo + $1.0)], 


and hence that 
| (1 — 02 + fF — fued] exp [—e( fb + o1)] de 


= {1 + O(f.'”’)] [ e(1 — [2 + f.F — fredo] exp [—€(fabo + ¢1,0)] de (51’) 


“0 








an 1 F 
= (1 + O(f.'”")] -—— - 
. ‘ )] Sabo Po 
One thus obtains the expression being sought, namely, 
aa 1 (3 + Yale)" l = fe) . a . 
(%) = y (1 — uz)"(1 — Bu.) \1 — y'7u/Sf, 1 — y'"u./fi Sabo ’ (52) 


where (¢) denotes the asymptotic expression for go . 
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It should be noted that the following is true: 


1/2 1/2 . 
vy “u,/ I Ua! f Bu, Qo rs 
hi = > — 172 ore... + | = a me ig = / ayes: 53 
: 1+ yu, i 1 — ¥'"u./f,— i. y(1 — w3)(1 — Bu,) (53) 
and consequently 
| — Bu, F (: + 7U, fa)" 1 — Bu, y" (59° 
(o,) = - : . ws ‘ 52 
sf Qo 3 L—Y¥y Ue/ Ji 1 — 7 “Ua/ Ji 


In order to check Eq. (52’), we compare the value of g» given in Sec. 10.1 for M, = 2, 


8 = 0.1,f, = 20, y = 1, with that obtained from Kq. (52’). 


One sees that 


0.4 (0.4)(0.1) 
== i = 0.374. , = a = 0.0605 
“1 ™ 2(0.66667)(1.1) ™ “2 ™ (0.66667) (0.99) aad 
and hence from the second of Eqs. (49) that 
1.66667 0.9333 
F = 0.273 log = + 0.0605 log = 0.501. 
£9 8 °0).3333 5 0.3333 


Remembering that 


1+ u,\*" ‘1.66667 \*:**** (3 = Bu. \**% (2:9388 ) — 
= - = 6494.3, = | —— = 3.4834, 
( ( a) setae j- & ) 0.3333 sane 


l— u, 0.3333 
it follows that the value of (go) given by Eq. (52’) is 
0.9333 0.501 


fy.) = ' :494.3)(3.4834) = 3290 instead of 268 
(%o) = 546 20 (6494.3) (3.4834) 3290 instead of 2681. 


One may improve the approximation by taking 0.99 u, instead of u, in the expo- 
nentials in the right-hand side of Eq. (52’); in fact, this assumption leads to the following 
value for (go) : (g) = 2820. Analogously, for 4, = 2, 6 = 0.1, f, = 23.6666, one obtains 
(go) = 14640 instead of g = 14707, etc. 

10.5. Justification of the logarithmic drag law. By means of Eq. (52’) one can easily 
justify the logarithmic law for the drag, as given in Sec. 10.1. Remembering the first of 


Eqs. (32) in Sec. 7, we may obtain from Eq. (30) of the same section that 


— si 0 > ool — Bu, 
dk, = —- S 2Qy°c, (1 — us) (go) de, = — = 2y°c, (1 — ty kf 
b : b Ao 
ec,” de (1 + yu. 2)" ( 1 — Bu, y" l+ncy moe oe 
7¢ i72 aa 2 ~_ i ° Team z 
(2y)'’"(1 — uz) 1 — yu I \l — y/ Ua/ Ji , ™ 
or 
dR, 2’ : 0 r 2\0 —n/ 7 
lt St we Hy] — ff)" — Bu 
de agb 


(54) 





Ly“ /7,\"" 1 —m, \" 
xX we 1/72. /¢ a eo . 
1 — ¥ Ua; Ji = J Ua/ fi 
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We find from Eq. (54) that 


es? = Lrg bog ( ct?” Se) — ng 
e270 — w]e ae, [2y1 — w)]'7F 
_ (55) 
9i/2 J per 
x log E yo" — way" aur | = A log (-2" aR) + B, 
ab dc, 


which has precisely the form of Eq. (38) with the added distinction that it now shows 
the influence that the Mach number and the thermic transmission through the wall 


have on the drag. 
We shall now compare the values of A and of B given by Eq. (55) with those obtained 


in Sec. 10.1: 


M, = 1.5; 6 =0.1 ; A = 1.835. instead of 1.781; B = 3.04 instead of 3.8 
=15; =0.25; A =1.74 “ “1671; B= 2.96 a agg 
M, =2 B=0.1; A= 1.89 “« “1,846; B = 2.92 = “3a 
=2: 2=0.25; A =1.73 “ 71698: B=264 “ “3.297 


Finally, it should be noticed that from Eq. (36) in Sec. 9 the following result may 


be derived: 


| log (c!"R,) lo E “0-8-8 

— > 73 og ae Se ee eee ——_ 
Oy — wr 8 §) ~— fod — wa) PF EL” 

(56) 


x (1 — a2)" “| = D log (Ryc}”) + E 


where 1 has been neglected in comparison with (b/f,)R; . 
On account of Eq. (42) in Sec. 10.3 and because of the expressions for A, B, D, E 
given by Eq. (55) and Eq. (56), it is consequently found that 


Ao ee » 
= 2y(1 on u;)(1 — Bu,)F *” . (57) 
which clearly points out the increase in thickness of the boundary layer effected by an 
increase in Mach number. 
HEAT TRANSMISSION BETWEEN FLUID AND PLATE 


The amount of heat which is transferred from the fluid to the plate, per unit length, 


is given by 
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a9 f° —l-, l= 7 l=, 
q=-5 | sofi+to-(i+)0)] ay 


(58) 
ie = — 
=—-—] sUB(U — U,) dy = Br, 
OX Jo 
where B is obtained from Eq. (10). 
Therefore, what has already been stated in Sec. 5 appears to be true, i.e., that if 
q = O then B = 0. - 
On the other hand (see Ka4rman’*) g = \,(dT/dy), , where \, denotes the coefficient 
of thermal conductivity at the wall; and since from Eq. (11) one has that (d7T/dy), = 


T,/U.)(1 — 1/y)(dU/dy), while (dU/dy), = r/up = paUcc,/2Qu, , it is found that 
fy Pp p 
T.¥—-1lie= ae eed. 
q=i\, => —— = — Us, = T.p.U. ~—— e, + - 
“U. v 2m, 2y My 


The total amount of heat transferred from the fluid to the wall, for a length x of the 
plate is then 

i Rem AVSY— 1s sen p Ya lz ne 

Q=] qdr=r, "75 —C, =,(1 — u2)"y'ToR: — C,, (59) 

Jo 2y 


which is consistent with the formula given by Karman in his paper,” except for the 
different expression for C, . 


APPENDIX 


11. Sample calculation of y, for MV, = 1.5, 8 = 0, and f, = 20. 


Jo (1 l 22 


= eer \l — uo Qu. -1— ue 
where 
Mevelcd” 
ul. = ree ; 
1+ 2/M-(k — 1) 
= 0.5571 
and 
— 7.1800. 


The numerical work necessary for evaluating the integral is carried out in Table IV. 
Column 17 gives the value of the integrand for the various values of the integration 
variable. Summing this column according to the trapezoidal rule gives 


go = 934.0. 
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TABLE IV. 
_ 2 3 | | 5 | 6 | 7 8 9 
— —— ——— — 
eles es @? fio | @? @/® |1+@/11-@]|! O/® 
0.0} 1.0000 0 | 0 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 
0.1 | 0.9000 | 0.0557 0.0031 | 0.9969 | 0.9938 | 0.9056 | 1.0557 | 0.9443 | 1.1180 
0.2 | 0.8000 | 0.1114 | 0.0124 | 0.9876 | 0.9754 | 0.8202; 1.1114] 0.8886 1.2507 
0.3 | 0.7000 | 0.1671 | 0.0279 | 0.9721 | 0.9450 | 0.7407 1.1671 0.8329 1.4012 
0.4 | 0.6000 | 0.2228} 0.0497 | 0.9503 | 0.9031 | 0.6644 | 1.2228 | 0.7772 | 1.5733 
0.5 | 0.5000 | 0.2786 | 0.0776 | 0.9224 | 0.8508 | 0.5877 | 1.2786 | 0.7214 | 1.7724 
0.6 | 0.4000 | 0.3343 | 0.1117) 0.8883 | 0.7891 | 0.5069 | 1.3343 | 0.6657 | 2.0055 
0.65) 0.3500 | 0.3621 | 0.1311 | 0.8689 | 0.7550 | 0.4636 | 1.3621 | 0.6379 | 2.1353 
0.70} 0.3000 | 0.3900 | 0.1521 | 0.8479 | 0.7189 | 0.4173 | 1.3900 | 0.6100 | 2.2787 
0.75} 0.2500 | 0.4178 | 0.1746 | 0.8254 | 0.6813 | 0.3669 | 1.4178 | 0.5822 | 2.4352 
0.80, 0.2000} 0.4457| 0.1987 | 0.8013} 0.6421 | 0.3115 | 1.4457 | 0.5543 | 2.6082 
0.85) 0.1500 | 0.4735 0.2243 | 0.7757 | 0.6017 | 0.2493 | 1.4735 | 0.5265 | 2.7987 
0.90) 0.1000 | 0.9014] 0.2514} 0.7486 | 0.5604 | 0.1784] 1.5014 | 0.4986 | 3.0112 
0.95} 0.0500 | 0.5292 | 0.2801 | 0.7199 | 0.5183 | 0.0965 | 1.5292 | 0.4708 | 3.2481 
1.00 0 0.5571 | 0.3104 | 0.6896 | 0.4755 | O | 1.5571 | 0.4429 | 3.5157 
| | | | | 
10 | 11 | 12 | 13 14 | 15 16 17 
o logio @ 7.18 40 logic! @ ® X ® o (3 2.3026 @ | 2+ @ @ xX @ 
0.0 0 0 | 1.0000} 1.0000} 0 | © | 2.0000 0 
0.1 | 0.04844 | 0.3478 | 2.2274 | 2.0171 0.2017 | 0.8008 | 2.8008 0.5649 
0.2 | 0.09716 | 0.6976 4.9842 | 4.0880 0.8176 | 1.6063 3.6063 2.948 
0.3 | 0.14650] 1.0519 | 11.2695 | 8.3473 | 2.5042] 2.4221 | 4.4221 11.07 
0.4 | 0.19681 | 1.4131 | 25.8882] 17.2001} 6.8800| 3.2538 | 5.2538 | 36.15 
0.5 | 0.24856 | 1.7847 60.9114 | 35.7976 | 17.8988 | 4.1095 6.1095 | 111.80 
0.6 | 0.30198 | 2.1682 | 147.3000 | 74.6664 | 44.7998 | 4.9925 | 6.9925 | 313.00 
0.65 | 0.32946 | 2.3655 | 232.0053 | 107.5577 | 69.9125 | 5.4468 | 7.4468 | 520.00 
0.70 | 0.35768 | 2.5681 369.9167 | 154.3662 | 108.0563 | 5.9133 7.9133 | 855.00 
0.75 | 0.38654 | 2.7754 | 596.2143 | 218.7510 | 164.0632 | 6.3906 8.3906 | 1376.00 
0.80 | 0.41634 2.9893 | 975.6600 | 303.9181 | 243.1345 | 6.8832 8.8832 | 2159.00 
0.85 | 0.44696 | 3.2092 1618.8148 | 403.5705 | 343.0349 | 7.3895 9.3895 | 2534.00 
0.90 | 0.47874 | 3.4374 | 2737.8125 | 488.4258 | 439.5832 | 7.9150 9.9150 | 4358.00 
0.95 | 0.51163 | 3.6735 | 4715.2222 | 455.0189 | 432.2680 | 8.4586 | 10.4586 | 4520.0 
1.00 | 0.54601 3.9204 | 8325.3333 | 0 | 0 | 9.0271 | 11.0271 0 
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THE SPECTRUM OF FREQUENCY-MODULATED WAVES 
AFTER RECEPTION IN RANDOM NOISE—II* 


BY 
DAVID MIDDLETON 


Cruft Laboratory, Harvard Universiiy 


Part I: Introduction and Discussion 


In our earlier paper [1]? expressions for the spectrum of an FM wave received in 
the presence of random (fluctuation) noise were obtained in the two extreme cases of 
no limiting and “super’’-limiting.** It is the purpose here to develop the completely 
general theory for the demodulated wave, taking into account the effects of arbitrary 
amounts of limiting on the noise and signal spectra and power of the low-frequency 
output of the receiver. We are interested in the spectrum because the signal-to-noise 
ratio determined at the output depends noticeably on the spectral distribution after 
demodulation when broad-band FM (in which the IF filter width is large compared 
with the audio response) is used. In the case of narrow-band FM (where the IF width 
is comparable with the audio) we are concerned mainly with the integrated spectrum or 
power output, so that spectral shape is not so significant (cf. [2]). The general theory 
is outlined in Part II, including an examination of important special cases, and a de- 
tailed study of noise alone is given in Part III; the Gaussian filter amplitude response 
is assumed. Remarks on signal and noise are contained in Part IV; applications to 
evaluation of the signal-to-noise ratio are discussed more fully in [2]. 

A satisfactory analytical model of the actual nonlinear elements in the receiver— 
the limiter followed by the discriminator—can be constructed if we assume: (1) that 
the physical discriminator is replaced by an “ideal’’ one which responds everywhere 
linearly with frequency; accordingly, the output current (or voltage) is directly pro- 
portional to the instantaneous difference frequency between the wave and the central 
or resonant frequency of the (symmetrical) IF, limiter, and discriminator bands; (2) 
that the filter response of the limiter is taken to be wide enough to pass the IF portion 
of the limited signal and noise without distortion due to frequency selection. In practise 
this means that a limiter band width several times the IF spread is needed; it can easily 
be obtained, since the limiter circuit is of necessity a low Q device. (If it were not, 
filtering would restore randomness to the noise, and the limiting would be nullified.) 
The case of the discriminator’s response is more critical, but if the linear portion of 
the actual characteristic is at least twice the r-m-s frequency deviation, distortion will 
not be serious, and our idealized model is then a satisfactory substitute (cf. [3, Chs. 4 
and 5] for a treatment of an actual discriminator when there is no limiting). 

*Received December 31, 1948. The research reported in this document was made possible through 
support extended Cruft Laboratory, Harvard University, jointly by the Navy Department (Office of 
Naval Research), the Signal Corps of the U. S. Army, and the U. 8. Air Force, under ONR contract 
N5-ori-76, T.O.I. The author wishes to thank Mrs. Roger Stokey and Miss Marilyn Lang of the Elec- 
tronics Research Laboratory, who performed the calculations for Figs. 1-6. 

tNumbers in brackets refer to the bibliography at the end of the paper. 

**In [1] is listed previous work on this problem, of which the papers by Stumpers [Proc. I.R.E. 36, 

1080 (1948)], Blachman [3], and Rice [9] are particularly to be noted. 
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The narrow-band wave leaving the IF and entering the limiter may be represented as 

V(t) = R(t) cos [wt + 0(t)], (1.1) 

where F is the envelope and 6 a phase angle. Now {(=w,/2z7) is the central frequency 

of the IF, limiter, and discriminator; both R and @ are slowly varying functions of the 

time compared with wot. Letting f(iz) be the Fourier transform of the limiter’s dynamic 
characteristic g(V), we may write the output of the limiter as (ef. [1)), 


Vo) = os I. f(iz) dz exp [izR cos wot + 8)] 


II 


> B,(R) cos [n(wot + 8)], (1.2) 
where 


n 
ae dt €&, 


BR) = 5 | f(iz)J,(Rz) de, (1. 
“C 


Oo 
— 


or 


(n = 0, 1, 2, ---). 


The contour C extends along the real axis from — ~ to + © and is indented downward 


about a possible singularity at z = 0. The various B,(R) are the envelopes of the n 
spectral zones produced in the limiter by its nonlinear action [4, Sec. 3], while the 
n@ (n > 0) are the respective phases. Only the band concentrated about f,(m = 1) is 


passed to enter the discriminator. One may show [1, Eqs. (1.7)-(1.9)] that the low- 
frequency output of the discriminator is 
. os. oF .. sila 
E,(t) = B,(R)é = f(tz) J (Rez) dz. (1.4) 
T JC 
The transform of the characteristic of our idealized limiter [1, Fig. 1] is 


f(iz) = 2e[1 — exp (—7R,z)]/(iz)’, (1.5) 


in which £ is a tube constant and R, is the level at which limiting takes place; the factor 
2 arises because both positive and negative portions of the wave contribute to the 
envelope. Explicitly, Z,(t) is found for our particular choice of f(z) [or g(V)] to be 


R, o<2 <R,, 
AR, = zea 

° —— ff a 2. 9: § /9- We a 

E,(t) = Bx04 x Fi(—1 1/2; 3/2; Ro/K’) (1.6) 
_ 2RJR, 
a T \R 





°(1 — Ro/R’)'? + sin™ Bel, R, < R. 
R 
The details of the analysis for the correlation function, the spectrum, and power follow 
in Parts II and III, and examples are illustrated in Figs. 1-6. 
At this point it is convenient to list the principal parameters that appear throughout 
the text and in the figures: b) = mean input noise power, at the IF output, Ay = peak 
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carrier amplitude, Ry = amplitude at which limiting takes place, R = envelope of the 
wave leaving the IF stage, p = Ao/2b) = ratio of mean carrier to mean noise power 
at IF output, tr = R,/(2b.)'”* = ratio of the r-m-s clipping level to the r-m-s noise 
level, w, = (angular) frequency displacement from exact tuning, w, = (angular) fre- 
quency proportional to the width of the IF spectrum (see Eq. (3.3a)), 2 = w/w, = 
f/f, = a normalized frequency, measured in terms of the IF spectral width, w) = 
maximum spectral intensity of the noise at the IF output, 8, x = limiter and discrimi- 
nator circuit constants, Do(t) = (angular) frequency modulation, R(t), r(t) = correlation 
functions, 7(¢) = phase of the modulated carrier wave. 

A number of general observations can be made. First, as shown in Figs. 1 and 2, 
the total power output drops with increased limiting, since less of the original wave is 
then passed. When there is no carrier but only fluctuation noise, we see from Fig. 3 a 
similar behavior, noting that in all cases when the limiting threshold (Ro) is very large, 
the output power remains constant for a given input noise voltage and varies directly 
with the amount of signal power, if a carrier is present. The signal is always suppressed 
when it is weak relative to the noise (cf. Eq. (2.10)). 

Second, as in any clipping or saturation process restricting the instantaneous ampli- 
tude of the disturbance, limiting spreads the spectrum; for noise alone, the precise 
extent of the spread is shown in Figs. 4 and 5. Figure 4 compares absolute values of 
the spectral intensity for varying amounts of limiting, and Fig. 5 gives a comparison 
of the same spectra, all now normalized to unity at f = 0, the point of maximum in- 
tensity. The reason for the spread lies in the fact that clipping produces an (infinitely) 
large number of new harmonics of the original wave, and the difference or beat fre- 
quencies between the original (adjacent) components, all relatively close to the resonant 
frequency fy , are the source of the added intensity in the low-frequency part of the 
spectrum. The effect, of course, is much more pronounced the heavier the clipping 
(R, — 0). For no limiting, the spectral intensity falls off — exp [—(f/f,)"], while in 
the instance of “super’’-clipping the decay goes as (f/f,)~*. The “tails” of the spectrum 
well away from zero frequency are quite extended (cf. Figs. 4 and 5 once more). One 
finds also, as earlier [1], that limiting, when the carrier is much greater than the noise, 
yields a vanishingly small noise spectrum at or near f = 0; if there is no limiting, how- 
ever, this is no longer true. Herein lies the superiority of broad-band FM over AM 
and narrow-band FM for large carrier amplitudes. The present paper, in conjunction 
with the results of [1] and [2], completes our theoretical discussion of FM reception in 


fluctuation noise. 


Part II: General Theory of Arbitrary Limiting 


As before (cf. [1]) the correlation function of the low-frequency output of the dis- 


criminator is* 





R(t) = [E(t)E(to + Day = «[B(R,)B,(R2) 6,8 2Jav « (2.1) 


The bar denotes the statistical average over the random variables describing the noise, 


*Throughout this paper we write the correlation function as Ro(t), and its normalized version as ro(t), 
to facilitate the distinction between the correlation function and the envelope R, and clipping level Ro . 
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Av 


while [ ], 


By a well-known theorem (cf. [5], [6]*) the mean power spectrum is the cosine Fourier 


indicates the average to be taken over the phases of the modulation, if any. 
transform of the correlation function: 


Wo(f) = 4 | R(t) cos wt dt, (w = 2rf), 
and conversely (2.2) 


R,(b) | W (Sf) COS wl df. 


The parrow-band input voltage V,(¢) of the limiter-discriminator element and the 
final low-frequency output /,(é) of the discriminator are represented precisely as in the 
earlier paper [1, Eqs. (2.3)-(2.11)]. The present analysis proceeds in the same fashion 
also, except that instead of the integrals A, and K, [l, Eqs. (2.20)] one encounters a 


more general version 


i y J i(é a+y : a : 
A | dx | dy a ‘ae exp (—irz — tyz2’), 
‘ ‘ ~ — 7 
(2.3 
; F F J (tx? L 4°)’ ; . 
K | dx | dy y ———— exp (—22z — iy2’), 
‘ ‘ o a ae 


which it is convenient to reduce to a single integral. Here é is a variable corresponding 
to the variable of integration in (1.4) which is introduced to account for the general 
dynamic characteristic [whose Fourier transform is f(7~)] considered here. As before, 
the following important quantities are needed: 


b, | (w — wo)"w(f) df; ?,(t) = = | w(f) cos (w — w,)t df; 
e tf /0 


(2.4) 


r,At) = (0 (0). 


Here w(f) is the mean input power spectrum of the noise, essentially determined by 
the IF filter response, which is assumed to be symmetrical about f, , the resonant fre- 
quency. The quantity bo is accordingly the mean input noise power. 

A more detailed discussion of the principal steps in the analysis is given in [1]. The 
final result is the low-frequency correlation function 


( 
— tk) €; rn 
R(t)» } 7 (f)° \ —¢,(b) > [3Ce.om.e+1 COS (kK + 1)(m2 — m) 
j Pages 4 \ 2 
+ 5Ci om.jz-1; COS (kK — 1)(m2 — 0;)] 
eee j Bo ; , 
+ oi(t) 5 | Heomsi,e COS k(n2 — m) — 5 We.om+i,e+2 COS (k + 2)(n2 — m:) 
er ans 
a Ky ome1.jk—-2| COS (kK — 2)(m2 — 2) 


*See also 8S. O. Rice, Bell System Tech. J. 23, 282 (1944) and [4] for further references. 
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(m + 1) —_ 
+ (2byp) ’*d,(t) mn = ss (da — Bla sn 0+1lhe.enst.b00 
— x, .) sin (k + 1)(m2 — m) 
= 5, es L 1 (3, Q2m+1,k ~~ Ki. om+i, I »;) sin | k _— 1 | (no — m)] 


¢ : E ce 2 I 
+ 2bopm n2 COS k(ne no] S82 ane + 9 (He. ome+1 = Sasa 


in which 
1.2 = wyl ot V; 2 and <7 = Wa + Gis = Wa + Dots ,2); 


t (2.5a) 
p = A,/2b, and ¥ = / D,(t) dt. 


The result (2.5) is formally identical with our earlier expression for Ro(t)y (ef. [1, Eq. 
(2.28)]) where only the cases of no limiting (R, —~©) and extreme limiting (Ry — 0) 
were considered. However, the amplitude functions 3 now take the more complex form 


; i (2 K ce) 
Rs snsn. dP ke) = = — —— [ f(ié) dé [ J ,(ré) dr 
= 20°"*“2 mim + 1)? ® Je" ~~ : 
: 2.6 
x | oP J(rp)S_(Aop) exp (—p*bo/2) do, _ 
(zu = 2m+n+k = 1,3, 5, ---). 
For the linear limiter considered specifically here, 3 becomes (cf. Appendix I) 
. sepia : xBp*’*(2/b,)"”” 
J _= — = ——£ —__ eee —~--= 
a Q'2m+)/2 mim + ky’? * qin'’*[m\(m + k)!]'" 
(2.7) 
1 ss’ r(—srd — sr’ "Te +s) , 
Oni | (s+ 1/2) + NTrB 2— 5 iFi(v +s8s;q + 1; — p) ds, 
where Tf = Ro/(2bo)'*, v = (u + q + 1)/2 = 1, 2, 3, --- , I, is given explicitly in 


(A1.15). [Similar results may be obtained in the same fashion for other types of limiter 
response, i.e. choices of f(zé) other than Eq. (1.5)]. 

As before, the complete correlation function R,(t) is obtained after the average over 
the phases of the modulation has been performed. The mean power spectrum follows 
at once from the Fourier transform of R,(t), according to (2.2). We observe that the 
output, as in amplitude modulation (cf. [7]), consists of three classes of modulation 
product: (1) (n X n) noise, produced by the beating of noise components with one 
another; (2) (s X n) noise, which is the result of signal and noise beating together; and 
(3) (s X s) signal harmonics. In (2.5) only the terms in the first bracket [ ] for which 
qg=\|k—11,k,|k — 2|, ete. = 0 represent (n X n) noise; the remaining terms of 
the first and all those in the second and third brackets [ ] are (s X m) noise, except 
in the latter when m = k = 0, in which case we have the signal components. 
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1. Signal Output. This follows immediately from (2.5) and is 
Ro(t)cexe) = Adlmne]avICo1 (2.8a) 
16 > 2 Dor . 9 l = >. 4 = 
= = 6 R3{ [ws + Do(t)][ws + Do(to + hap — 2D) 
- 
-- ’ (-2) (—r5)' = ([ vo + ¥(l + 1)/2 — ¥(8/2 — I)/2 — logr, (2.8b) 
er 2/, 2i+ 1 , ao 
1 1 (ll + 1),¥(l +n + 1)(—p)"\’ 
ee FS ee ee Se SET, Mie A dl AF a Sone 
2i(21 + >| STE ws a (2),n! 
10 10gig A 
(0) A= Ro (t)isxs)/<1 Nerv. 
(b) = Ro(O)isxsy/<(7Q)* > gy, (POWER) 
(c) = Ro (O)isxsy/ <(1)? oy. = (0.C.)? 
db 
p= ak/2bo 
% = Ro/, 2bdo 
-10 
-20 
10 logig P 
aa “10 0.0 10.0 20.0 300 
Fic. 1. Discriminator output—(a) Correlation function for signal components, (b) Mean signal power 


output, (c) Mean square d. c. output—as a function of carrier strength for various degrees of limiting. 


Here the abscissa is 10 logio p, where p = (mean carrier power from JF)/(mean noise power from IF); 


7 is the modulation, including any deviation (wg) from exact tuning (fo = f-). 


from (2.5a) and (A1.15); Ww is the logarithmic derivative of the Gamma-function. Note 
that because of our assumption of an idealized, linear discriminator, the modulation is 
received undistorted. Figure 1 illustrates Ro(t),.x.) as a function of p for various de- 
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grees of limiting. When there is no limiting at all (r, ©), or when the limiting is ex- 
treme (r, — 0), we obtain (cf. Eq. (A1.18)) the expressions of our earlier analysis 
[1, Eqs. (2.36), (2.37)]. In general, when the carrier is very strong compared with the 
noise (p ©), we find (cf. (2.8a) and (2.19)) that for all degrees of limiting such that 
R, < Ay (orto < p’””), 


16 


x 


Rolf) cex2)]p+0 «BE ni nelav R; 2F'(—1/2, 1/2; 3/2; Rj/ Ao)” 


r 2 
4 f6'TivisleR'([1 — Rig”? + A%sin-* Be), sea 
T Ry Ay 


| < Ay . 


Here the noise is suppressed by the signal. Compare (2.9) with the square of (1.6); 
the dependence on R, and A, (=R if the noise is negligible) is the same, as we would 


expect. 
At the other extreme of carriers weak compared with the noise (p — 0), (2.8) re- 


duces to 


(1), o)-o te (4 gam.) forge).p” (-UBf—B" ( 2m 
Ro(®)¢ex2]y-0 = (! 8aR,) Civile _ > ni? a ag 


m=1 





(2.10) 
x (v(m) — ¥(38/2 — m)/2 — logr, — 1/2m(2m + oi} , 


This shows at once how the stronger noise suppresses the signal: instead of the (power) 
output being proportional to p, it is proportional to p’, just as in the analogous situation 
for the detection of AM signals by a half-wave rectifier. 

2. The d-c output. The d-c output is obtained from the correlation function if we let 
t —-~ in (2.5). The square root of the result gives us the desired mean amplitude, 


which is 
(Eo(t)Jav = KBAolwa +- Do(t)JavRoor . (2.11) 


By a different method (cf. [2]), we may express the d-c as an integral, 


a@ 


[EvOlee = 2unp'e? | By(r)I,(2rp) exp (—r*) dr, (2.12) 


where B,(r) is given by (1.3) on replacing R by r and multiplying by (2bo)'’*. Numerical 
results are also shown in Fig. 1 (for the square of the d.c.), which may be obtained 
either by summing the series in (2.8b) or by graphical methods applied to (2.12). Simple 
forms of [£,(2)],, occur when there is either no limiting or extreme clipping, (cf. [1, 
Eqs. (2.36), (2.37)]). 

3. Mean power output. By a well-known theorem [4, Sec. 2 and Appendix II] the 
mean power may be derived from the correlation function on letting t — 0 (ef. (2.2)). 
The result in this instance is the double series given in [1, Eq. (2.33)], on replacing H 
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by &. Ag: 
namely, 


Wo 


The curves 


levels in th 
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[2]), 


1in, a simpler expression has been found in the form of an integral (cf 


no 


K exp ( py pli | B,(r)* exp (—r*)[Io(2rp'”*) + 1,(2rp'”’)] dr/r 
(2.13) 
¢(0) f° » - . 
7 ! | B,(r)° exp (—1°)I,(2rp'””) dr/r 
yo Jo 
of Fig. 2 show how W, “tide nds on carrier strength for various clipping 


e important special cases [n lav = O (no modulation), and [n” 7. —,(0)/2bo . 



































































































en. ee | “of ae « bf oF a T 
[ | | | ; | | | 
ac i a a we | 
ae Py, [KN2>=be/2 
9:0 [ Wo /(K®42b,)] —F t t ot 7 ‘ - do] 
20 eS {+ _+ _} Bez } | | 4 ; 
| p=As/2b, ) a SGE25 .2 a eae 
b- 1 = Ro/V20o - } am 
| | | | | 
ERGhAEERSET AER aan 
—— MODULATION | | 
—-- NO MODULATION 
10 f 
db 
° 
- 
20 r - 
) - +++ 
| To,<TE2= 0 |e? a 4 S| } NI 
+ + + + + + ——— N + + t + + _ — I< t 
| | } | i" », | \ | | ‘ | NY 
+—t {——+ —4—— 4 4 } + = 
, 4 Sega et | N{ ae } JN W ~~ 
+—+ . + i 4 + +\ + \ } 7 WHO’ +__ 
rT rrr | ~ wRROSS 
-30 =: 4 4 + 4 4 4 4 + a ae —_j—_»,—++ SS = ~j— 
| | | | | } | ‘i be \ | \ 
COT saa Beaceeeach 
| | \ | 
-——+ + t — +4-- <a | —-+——- + --—- | _ T og ; a TR po — 
| | } | \ \ | SJ 
+ t - — —_ —— —+ T + T — ko $+. 
| | 10 logig P | =| db mk i | 
a" + ~¢—~ 4 4 —4-—~-4--- -— ay 4——_+4 
| | ff | | | | retry PT SX 
-40 | | = | i coal ae | 
20 “10 r) 10 20 30 
Fic. 2. Mean power output of discriminator as a function of carrier strength, with different degrees of 


limiting. 
obtained fror 


Her 


abs and 7 is the modulation, if any. The quantity b» 


) 


issa is 10 logio p (cf. Fig. 1 
, represents the second moment of the input noise spectrum w(/) 


1 
e the 


n Eq. { 
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Specific calculations were made with the help of Simpson’s rule (ef. [2, Part II]). In a 
few special cases, however, W, assumes a closed form 


(no limiting: r, ~«©)W, = 2byx's| alia 0 ee =P) a #00) | (2.14) 
Io 


2p 
and 
° - » o 4xB E *9 J 
(extreme limiting: r, — 0)W, = —f,] bo exp (—p)\ [pln lew > $2(0) /bo] 
vg 
xX [y — Ei(—1o) 
(2.15) 

+- 2 n’ le p”'*(2m + 3) — o.(0)p"*'(m + 2)/bo]/(m + 1)(m + 2)3}, 

where y = 0.2731 and Ei (—r;) is the exponential integral 


a, 
_ le e” dy/y. 
When the carrier is very strong we obtain (2.9), (¢ = 0). (Other limiting forms of the 
above are given in [2, Part II]). 

4. Two important limiting cases. 

A. Strong carrier with limiting (Ag ~©; Ry < Ao , bo # 0). Here the noise and 
limiting levels are assumed small compared with the peak carrier amplitude, i.e. 
Ay >> (2b))'*, Ayo > Ry [however, it is not necessary to assume that R, > (2b))'”’). 
Then, in our expression for the output correlation function, obtained from (2.5) we 
retain only terms in 6? and by , and discard the higher powers of r>' and p~'. Our correla- 


tion function reduces to (p © in what follows) 


R,(t)y 2 —H iP2(t) COS (n2 — m) + 4 Ag 5Coords (Et) -(m: T no) 


(Hore — Horo) sin (n2 — m) (2.16) 


2 


perce ae 
‘do(t) m2 COS (Ne = m) * (Hiroe ioe 5100) + Aomn23Coor ° 


* 2h 


Now we need the limiting form of 3C when p ©; this is most easily found from (2.7) 
and the asymptotic development of ,F, [4, Appendix IIT] to be 


F ‘ 8 (2) ‘ ” SBuaaitaicie 
J — o — 0 } 
ieee: T bo, [m!(m + k)!]'” 
(2.17) 
/ ' —sI(—s)*(ro/p)'T + s[1 + v+s)\v +s — g/pt+-::') de 
Qri (s + 1/2)T(s + 13/2 — s)l(¢q+1—v-—s8) ; 
Keeping only the leading term in the series, we easily find from (2.17) that 
2 2b,)'”” 
Hore —, Horo os 3oo01 ’ and IC 102 — IC 100 = (2bo) Hoo ’ (2.18) 


as Ao 
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where 
a KBR, 1 [f° T(—s)(Ro/A)” 
I, re | ds 


Aomr'”* Qri J_«; (8 + 1/2)1(3/2 — 8) ”’ 


_ 1xBRy oF, = 1 2. 1 2: 3/2; R3), (2.19) 
wAy 
9 r 
_ = |B (1 — R6/Ao)'? + sin” Re | (Ro < A)). 
re Ao Ay 


The general procedure for evaluating the contour integral for 3y9,(p — ©) is available 
in Appendix I. For the noise part of the output, the correlation function (2.16) is, there- 


fore, 
Ro(t)noiae 2 [—2(t) cos (m2 — :) + o,(8)(m + 2) sin (m2 — m) 
(2.20) 
+ o(t)min2 COs (ne — m))Hoo , (p—@; A, > R,). 
We observe that the coefficient of 3C*y9, in (2.20) is precisely that given in [1, Eq. (2.40)], 
when A = 2, i.e. when there is extreme limiting. Our more general result also includes 


the cases of moderate to negligible amounts of limiting (Ry — Ay , Ag > (2b,)'””). 
In a similar fashion (ef. [1, Eqs. (2.41)-(2.47)]) we obtain the alternative representa- 


' 0 0 0 ra] 
0 oise = 7 "001 eit o o - ‘0 
Folds : ‘( dl + 2\3, si os (Y 


x | cos E + f? Dt) ww | ae/xsp, 


J¢ 


tion 


in which 2o(=27/w,) is the period of the modulation and 0"¢,(t)/dt" = @,(t), (iq. (2.4)). 
Expansion of the integrand in a (cosine) Fourier series, using the integral form of @¢,(¢) 
and the Fourier transform relations between the correlation function and the mean 


power spectrum of the noise, gives us finally 


W(f)< (4x8) oF',(—1/2, 1/2; 3/2; Ro/ As) 


wA,y 
(2.22) 
> A,w™[wlwo + nos + ws + w) + ww. + nw, + w, — w)], 
n=0 
(Ry < Ao) 
Here A,, is the amplitude of the nth term in the development of [cos (yz — )],. 
Equation (2.22) shows that when the carrier is strong and there is limiting at or below 
the maximum carrier level, the noise spectrum always vanishes at f = 0 and is small 


in the vicinity of zero frequency. This accounts for the great improvement in the signal 
vs. the noise (ef. [2]) when broad-band FM (maximum modulation frequency small 
compared with the maximum carrier deviation) is used at large carrier levels with 


sufficient limiting. 
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The mean power output, including the signal, is found at once from (2.9) and (2.20) 
by setting ¢ 0. We have 
(p > 0)W, = R,(0) = {Az + ¢,(0)[n"].. — $2(0)} 520 
(2.23) 
= Woon {pla lav + 3([7 lav — 72(0))}, (Ry < Ap) 


showing how the strong carrier suppresses the noise [second term of (2.23)], as in the 
detection of amplitude-modulated waves. 

B. Strong carrier, little or no limiting (Ay -~&©, Ag < Roy , bo ¥ 0). In the extreme 
of exceedingly strong carriers, when A, < R, the situation corresponds in the first 
approximation to the case of essentially no limiting, for which 5C becomes from (A1.18), 


KB (2)? pS" r([2m + k +n + g)/2) 
nS is l(g—2m—k—-n+2/2) ° 4) 


The correlation function (2.16) then reduces to 
Ro(t)y = « B’{—¢.(t) cos (nx — m) + 2bopm ne}, (Ry, > Ap) (2.25) 


and the noise spectrum associated with the low-frequency output is again precisely our 
earlier result [1, Eq. (2.42)], viz: 


Wo(f) = «’B’ ) A,|[(nw, + ws + w) wlwy, + nw, + wy + w) 
aie (2.26) 


+ (nw, + wy — w)’wlw + nw, + w, — w)], 


where we note, in contrast with (2.22), that the spectral intensity at and near zero 
frequency does not vanish. Broad-band FM under these conditions is accordingly much 
less satisfactory (cf. [2]) than when limiting is used. The mean power may be obtained 


from (2.14) on letting p—-@. 


Part III: Noise Alone 


This is the simplest case to handle analytically; a complete discussion follows. Since 
there is no carrier, p = » = 7 = 0. We may therefore, from (2.5), write the correlation 
function of the low-frequency output of our discriminator as 


; = 9 2m+ ] ~ 2 2 m 
Rol t) — ; {2h ¢.(t)r,(t)” ' + 2 Io, 2m4 1,0 ¢,(t)’ro(t)’ 


(3.1) 
1 9 2 2m+2 
ie 9 Ko 2m +#1,0 o,(t) ro(t)” \ 


Again, in the special instances of no limiting (tr) ©) and super-limiting (t —0) (3.1) 
reduces to the much simpler forms of our earlier treatment.* A general expression for 
the mean output power W, is found by setting ¢ = 0 in (3.1) above; however, a compact 
form, from which calculations are more easily made, may be derived from (2.13) on 
letting p — 0, viz., 


*Cf. Eqs. (3.2) and (3.8), and Sec. 3(a) of [1] for details 
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| | a 
| B,(r)* exp (—*) ©, (3.2) 
= 


0 


p (0) 


Wolr-0 ‘i 


Figure 3 illustrates the variation of the mean noise power as the limiting threshold is 


changed. 
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Fia. 3. Mean low-frequency noise power output of the discriminator for different degrees of limiting 


when there is no carrier; b, is the second moment of the spectrum w(f) Eq. (2.4 


To obtain the spectrum we need only take the Fourier transform of (3.1) according 
to (2.2). In the specific calculations, we assume a Gaussian spectral distribution for the 
input noise, corresponding to the composite Il-limiter-discriminator frequency response. 


We have 
> ? 1/2 ‘ ‘ 
w(f) Wy exp [—(w Wy) /wy], and Wo Qhow ’“/wy . (3.3a) 


From (2.4) we find also that 


od, (t) h, exp (—w,l 1), ¢d, (0) Sw, t do(t) 
(3.3b) 
o (7) bs (wt Y. \p | za. r Tol: wr, az. 


The spectrum becomes finally 


f f ‘om 
W,(f) Or “bow - a: Jf (1 — /)(m + 2) ‘een + 2) 
m 
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f,, . Figures 4 and 5 show the spectra for different degrees of limiting. 
Figure 5 is the same as Fig. 4, but now the spectra are normalized in such a way as to 
have the same intensity at f 0. This demonstrates more effectively how limiting 


where 2 = w/w, = f 
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Fic. 4. The spectrum of the noise output of the discriminator for various degrees of limiting when 
there is no carrier (p = 0). 


spreads the spectrum. A table of spectral intensities and a brief outline of the method 


of calculation involved in Figs. 4 and 5 are given in Appendix II. 
With no clipping or, conversely, with very heavy clipping, the spectra are given by 


[1, Eqs. (3.7) and (3.9)]. These distributions are included in Figs. 4 and 5. It is interesting 
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to observe that well out on the ‘‘tail” of the spectrum, where 2 ~~, the intensity falls 
off in the following ways:* 
Wo(f)] 2 or KB bow, exp (— 2") 


, 


. | ee 
W(f)lr,-0 = — R2x°6%w,27'. 
Q-~ T 


Wo(f)/W,(0) 
| 


N = NORMALIZATION FACTOR 


Ty = Ro//2do 


———— INPUT SPECTRUM: e-.* 


1.0 “ee _ 
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FACTOR (N) #10 109i N 
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Fic. 5. Normalized spectrum of the noise output of the discriminator for various degrees of limiting when 
there is no carrier (p = 0). 


For intermediate amounts of limiting the behavior varies between the rapid falling off 
[— exp (—2°)] for no limiting when 0 is sufficiently great and the relatively slow decay 
(—Q"") of super-limiting. Note also that for large clipping levels (ry, ©) the spectrum 


*See Appendix II. 
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is independent of limiting, while for very intense clipping, the spectral intensity is 
proportional to rj (tr) — 0). 


Part IV: Remarks on Carriers and Noise 


When a carrier accompanies the noise in the receiver, our final low-frequency output 
is considerably modified, particularly if the carrier is intense relative to the noise. 
Modulation further distorts the continuous part of the output spectrum, so that a 
precise calculation of spectral shape becomes formidable indeed. In the more general 
cases including a signal we use (2.5) and average over the phases of the modulation, 
according to 

1 To 

R,(?) baa T, i Ato moaylro(t, bo)v = [Ro(Owlav 5 (4.1) 
to obtain the complete correlation function. Formally, we have only to replace the 
amplitude functions H of our earlier analysis [1] by the more general 3 (ef. (2.6), (2.7)), 
yielding the behavior for all degrees of limiting, the most important cases of which 
(r,) >, ry —> 0) have been discussed in the earlier paper [1, Secs. 3,4]. In the general 
instance of arbitrary limiting the results of [1, Figs. 5-11] are expected to apply quali- 
tatively, except that as the degree of limiting decreases, there is less spreading of the 
spectrum for a fixed carrier power (p = constant). Out on the tails of the spectrum 
(Q —) the intensity [for an originally Gaussian distribution of the type (3.3a)] falls 
off at least as exp (—’) for no limiting and at least as Q"’ when there is super-limiting. 
The presence of a carrier (p > 0) increases the rate at which the intensity diminishes, 
especially if the carrier is strong: the noise is then largely suppressed and the spectrum 
is proportional to 2’ exp (—”) (ef. [1, Fig. 8]). On the other hand, for weak carriers 
the noise in turn suppresses the signal (Eq. (2.10)), and the distribution of the spectrum 
will follow as if noise only were present (cf. Part III and Figs. 4 and 5). The effects of 
modulation do not become noticeable, of course, until we attain large carrier powers; 
Figs. 7-11 of [1] show typical distributions in the extreme of strong carriers with simple 
forms of modulation (ry —© or fm — 0). 


APPENDIX I: The Amplitude Functions 3C, on+.,.(P; Ro) 


The quantities #, which appear in (2.6) for the correlation function, become for 


the linear limiter used here 
(2m+k)/2 
_ 





RK = Hz, 2m+n.a(P; Ro) i 9 ’ (Al.1) 


ae 


(m+!) 20 mim + k) i’? 


where 





_ —2ixB [ (3 = exp, {—iRe) de [ J (rz) J (rp) dr 
x Jc z si 


(A1.2) 
x [ p**' J (Aop) exp (—bop’/2) dp, (a = 2m + k +n) 


in which, from (2.5), we note that u = 1, 3, 5, 7, -- 
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We start by integrating over r first, observing that this integral is a special case of 

Weber and Schafheitlin’s more general expression. Using Gegenbauer’s result* we obtain 
e Zp 1/6 es 23 2 272 ¢ 

L,(z, p) = | J (rz) J (rp) dr = = naa ol”,(3/4, 5/4; 2; 42’p°/[2” + p')’) (A1.3) 
2(z ae ) 
which is defined for all z # p. When z = p, 2F, is logarithmically divergent, and /, 
assumes the formal value I'(0)/2zz. Since the singularity is only logarithmic, it is safely 
removed in the following integration over z. We have** 


“0 


L,(z, p) dz 


f 1 — exp (—Mee 
La(p) [ l — exp. tRoZ) 
c 

(A1.4) 


= 4, $F B/NA6/9).2" f° we sin Ropu 4, 
2 (2),,n! Je (xy? + | ial 


Next, we represent sin R,pu as a contour integral involving I-functions, reverse the 

order of integration, and finally sum the series. Thus we write 

1/2 m4 2s+1 , 
: a R, pu I'(—s) " 
sin R,pu = —— [ ( : ) ——_—-— qjg, (Rop > 0), (A1.5) 
2nt Jo: 2 I'(s + 3/2) 
where the contour extends along the imaginary axis and includes the pole at the origin. 
With the aid of Cauchy’s theorem the integral may be evaluated to give sin Ropu, by 
swinging the contour around in an infinite arc in the half of the complex plane for which 
R(s) > 0, and hence enclosing all the singularities of ['(—s). These singularities are 
simple poles at s = m (m = 0, 1, 2, ---), whose residues are (—1)”"/m! The convergence 
of the integral and the proper vanishing of the integrand as | s| —© along the infinite 
are are easily shown with the help of the asymptotic expansion of the I'-function,f viz., 
l(+z+ a)“ exp {(4+z2 +a — 1/2) log (+z) z + (log 2x) /2 + O(2z)}, (A1.6) 
| arg (+z + a)| < z, | arg (+z)| < a, O(z) > 0, |z| +o, 


For the integral in (A1.4) we find that 


[ sin R,pu in 
ne du 
7 I+u)”" 
1/2 hk | 28+1 + s 8 aco 2n+2s8 
.— | (Bee) F I ~ ds | = —— qdy, (A1.7) 
mat Jaws \ @ P(s + 3/2) Jo (x? 7 
dali [ (Bee)""" (s\n +s + 1/2)T(m —s + 1) ] 
= — . . - = sruer-ar -—agueeememmmt ( 
20(2n + 3/2) J_o; \ 2 I'(s + 3/2) 





*G. N. Watson, Theory of Bessel functions, 2nd ed., Cambridge University Press, 1945. See p. 407, 
Eq. (1), which is incorrectly given; a factor T'(A/2 + 1/2) should be inserted in the denominator. 
**Termwise integration is allowed, since the series is uniformly convergent, except at z = p, where, 


however, 
| | g(z) | dz exists. 
(z=p) 


See Titchmarch, Theory of functions, 2nd ed., Oxford University Press, 1939, p. 42. 
tSee E. T. Whittaker and G. N. Watson, Modern analysis, Cambridge University Press, 1940, p. 279. 


1950] SPECTRUM OF FREQUENCY-MODULATED WAVES 75 


on evaluating the Beta-function (which is defined, since R(s) < 2n + 1/2). Inasmuch 
as T(n + s + 1/2) (n — s+ 1) = I'(s + 1/2)T(1 — 8)(s + 1/2),(1 — 8), , and since 
(3/4),(5/4),2°""-' = [(2n + 3/2)/x'”, from the duplication formula for the Gamma- 
function, we see that the series in (A1.4) becomes 

1/2 


(s+ 1/2),(1 — 8), oy . athe J 
} ® (2) nl = Fis+1/2,1-8 2) =7oyppeery Als 





n=0 


by the well-known summation formula for ,F, . The integral L.(p) reduces finally to 


im? ff —sI'(—s)"(Rop/2)"*' ! 
(oe) = = 45— na sp. Al. 
he) = tat « @+1/2Te + DT@2—9~ ining 





Integration over p, with the help of [4, Eq. (A3.7)], gives us I, : 


pT (Aop)La(p) exp (— p bo /2) dp, 


Pi — 27 xB [ 


T 





KB p”~ & yn 1 [ ; —sI(—s)*r,"*'T + 8) 
! 2ri J_; (8 + 1/2)T(s + 1)T(3/2 — s) 


Tv q- bo 


(A1.10) 

xX ,Fi(v + 83 ¢ + 1; — p) ds, * v=(ut+ q + 1)/2 = 1, 2, -->), 
where p = A;/2b,) and ry = R,/(2b.)'”’. This form is particularly useful for large values 
of p, as then we may replace ,F, by its asymptotic development before evaluating the 


contour integral. 
We observe that the integrand of (A1.10) contains a simple pole at s = 0 and double 


poles at s = l(l = 1, 2, ---). We need, therefore, to evaluate an expression of the form 
1 ee | : ‘ 
L=5-—| 1(—s)*f(s) ds, (A111) 
an 


where f(s) contains no poles in the right complex half-plane [R(s) > 0] and I'(—s)f(s) 
is finite at s = 0. The residue R, at the double poles s = l(l > 1) of I'(—s)’ is found 
from the Laurent expansion to be 


R, = | ((—s + p*r(—9*s00) | (A1.12) 


s=I>1 
Using Euler’s representation of the I'-function, 


I'(—s) = lim [n!n“*/(—s)(—s + 1) --- (-s + n)], 


we get 
R, = lim ({—2n!’n-*' f(D) log n + nl’n7”' f(D) 


+ QnPn- f(D[—1/l + 1/(-1 + 1) 
(A1.13) 


fees $1/-14141/24--- +1/a— Dj} 


+ ((—9(—1 + 1) --- (—1)°-1°- --- @ — 9"). 
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J 
Lr) 


Since 
lim (1 aianiy see eo log n) 
Bae 2 3 n 


= y = 0.5772 --- , and ¥(m) = | +5t+a+ vee ti 


m 
where y is the Euler-Mascheroni cqnstant and y is the logarithmic derivative of the 
T-function, (A1l.11) becomes 


L = TO)fO) + > a [/'O — 24D f(). (Al.14) 


[The logarithmic derivative of the I'-function is tabulated in Jahnke and Emde, Tables 
of functions, Stechert & Co., 1938, pp. 9-22; note, however, that Jahnke and Emde’s 
definition is related to ours above by yy.n. (z) = ¥(z + 1)]. The leading term of (A1.14) 





represents the residue at s = 0. If we observe that 
d _ . —y(l + 1)  _ 4 W(3/2 — D) 
— Ts a ‘le ad 3/2—s pe oo Pe 
(a1 ee a Ta+1)?’ last 8/2—-9 f = T3/2- ) 


in ee 


{a} = 2(log ro)ro* 
and that 
FO Do + L+ n)(=p)" 





(3/2 — 1) = (—1)'n'”/2(—-1/2), , 


we obtain finally for (A1.10) 


: 4 op” 2 B/2 : 
r, = eT . I'v), Fv; g + 1; — p) 


)( = 1/2) (40) - (2) 
“—* (i) 21+ 1 


l= 


x ([ veo + Wl + 1)/2 — ¥@/2 — 0/2 - sp tas ins 


— log Fie + i¢q+1;—p) 


1s @+Duy vo + 1+ apy )} 
p> (¢ + 1),n! ’ 
v=mt+(kt+n+ q+ 1)/2 = 1, 2,3, --:). 


In cases of extreme or super-limiting (A1.15) immediately gives the behavior as a 
function of the relative limiting level r, as r,) — 0; we have accordingly, 
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4r5 “a(? y" pr((2m +k +n+ q+ 1)/2) 


do q\{m!(m + k) i}? 


lim 3Cz,2m+n,e(P} To) 
r.-0 Tv 


xX Fi((2m +k +n+9¢4+ 1)/2;¢+1; -pD 
(A1.16) 


4 
= Ro BKHy om+n.a(P} 2), 


where Hy.omsn.¢ is the amplitude function defined in an earlier paper [1]. At the other 
extreme of no limiting (r, —«) we attempt to expand (A1.15) or (A1.10) in a series 
in (tf) *. This is done by displacing the contour in (A1.10) to the left of the imaginary 
axis, a process that is valid for all R(s) < 0, since it is readily shown from (A1.6) that 
the integrand vanishes as exp (—7 | I,,(s) |/2), when J,,(s) —> +. We assume for the 
moment that v is not an integer. Then the first pole for R(s) < 0 occurs at s = —1/2, 
and succeeding poles lie ats = —v — l, (l = 0, 1, 2, ---). Applying Cauchy’s theorem 
to the region bounded by the displaced contour, the infinite arcs at I,,(s) ~ +, and 
the original contour, we obtain the asymptotic development 


— «Bp 2)" mo 27 Pe a : oe 
rs a! b, 9 Te 1/2)-,Fi@ — 1/2;q¢ + 1; — p) 
(A1.17) 

-1-TO)T? + Dr” < (v)3(v + 1).(t0)~" } 

+ nr °T(3/2+»)rd — v») w+ 1 — 1/2)(38/2 + »),l! Fk—be + 1; — wy. 
Now letting vy — 1, 2, --- we see that the series in (A1.17) vanishes (since ['(1 — v) > 
+o, or more fundamentally, as a consequence of ['(v + s)/T'(1 + s) not having poles 

when R(s) < 0), and we are left with 








xB 2 n/2_@/2 
lim Ky, 2m+n.c(P3 To) = 2 (2) a T((2m + k + n+ q]/2) 


Io @ 


(A1.18) 
X F,((2m + k + + gl/2;q +1; — p) = BeHi2men.c(D; 1), 
which agrees with Eq. (2.29) of our earlier work [1]. (If we note that (sin Rppu/u)|Ro.. = 
r5(u — 0) we obtain (A1.18) directly from (A1.4) and the integration over p, since 
only the term n = 0 yields a nonvanishing result). Unfortunately, the above approach 
does not give the succeeding terms in rj’, and other methods must be sought. 


APPENDIX II: Calculation for the Spectrum of Noise Alone 
The general spectrum of the noise output of our discriminator is given by (3.4) and 
the pertinent amplitude function 3 follow from (A1.15) when p = 0, namely, H,,2m,0 ; 


5o.2m+1.0 » IC2,2m+1,0 - Explicitly, these quantities are 


- 4“ 2 1/2 
IC; 2m,0 _ (m + 1) “Gn (To); Ko 2m+1,0 = (2) G,,(to); 
0 


1 1/2 2 1/2 
Ke, 2m+1,0 _ Gat (2. Gavi Ge), 


(A2.1) 
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where 


— 4 xBr, =. (—15)'(—1/2),(m + 1), ( 2! ) 
G(f>) = ; a 
m\Fo) rr {1+ p> (7!)' 21+ 1 


[y(l) + yl + 1)/2 — y(3/2 — 1)/2 — Wm + 14+ 1)/2 (A2.2) 


. log ro] 
. — = log ro]. 
2(2i +1)" 2°" 
The series for G,,(%) converges with satisfactory rapidity for all r, equal to or less than 
unity, but unfortunately the series (3.4) for the spectrum does not converge at all 
speedily when rp — 0.3 or more. Since accuracy to within a few per cent is quite sufficient 


{Wo (#)/ Wo(0)} 
A 
1>#R,//eb,; N= t/fo 
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The normalized spectrum of the noise output of the discriminator as a function of limiting level 


when there is no input carrier (p = 0). 


for all practical purposes, it is possible to get around the convergence difficulty by 

interpolation. To do so we first calculate data for Wo(f) in the interval 0.01 < ry < 0.3, 

which can be done precisely. Next, we observe that for values of ry — 2.5 or more the 

effect of limiting on the purely noise wave is slight. The probability that noise peaks 
—-6 


= ¢ = 


2.5 times in excess of the r-m-s value will occur is proportional to exp (—r) 
0.002, or about 0.2 per cent of the time—a negligible effect. The resulting spread in 
the noise spectrum due to clipping is correspondingly trivial (cf. [8]). Accordingly, we 
may safely use the values of the spectral intensity, when r, © [1, Eq. (3.7), when 
tr) > 2.5] leaving an interval of about one cycle over which to make the interpolation. 
Because we wish eventually to be able to normalize the spectra so that the maximum 
intensity is unity, and since Wo(f)... occurs at zero frequency when there is no carrier, 
we need make only one interpolation, i.e. for W,(0) over the range 0.3 < rm < 3.0. 


i et 
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From this data (correct to within two per cent) we normalize the spectrum at the calcu- 
lated points, obtaining the curves of Fig. 6. A short table of normalization factors N is 


given below: 




















TasB_e I. 
To | N* Ty N | | To N 
————ee {= 
| 
01 2.644-10-4 0.1 4.78-10°-? 0. 0.370 
.03 2.47-10-3 0.2 0.176 0.9 0.375 
05 [oe | 0.3 0.257 | 1.0 0.379 
.06 1.12-107? 0.4 0.303 1.5 0.382 
07 iw |< 0.5 0.332 2.0 0.386 
08 2.51 * 0.6 0.350 3.0 0.389 
.09 3.48 “ 0.7 0.362 ro) 0.3963 




















(*The values in bold face type are precise; N is the number by which Wo(f)norm. must be multiplied 


to give Wo(f absolute.) 

The limiting behavior of the spectrum at large frequencies is of considerable interest. 
It is possible to obtain at least the leading term in the asymptotic expansion quite 
simply by the technique of Appendix I (cf. Eqs. (A1.16)-(A1.18)). We begin first with 
the simple case of super-clipping, in which we have to examine the series (3.8) of [1]. 
We transform the series as follows: 


a exp [— 2°/2(m + 1)] _ > (— 92/2)" ‘ } 1 


m= (m + 1)°” rer ! mao (m + 1)**°? 


(A2.3) 
= D $(3/2 + k)(—9°/2)*/k!, 


k=0 . 
where ¢ is Riemann’s Zeta-function.* Next, we use a contour integral (cf. (A1.5)) in- 
volving the Gamma-function to represent the transformed series, viz: 
sine 


' I 
(3/2 + k)(—9°/2)*/k! = — 
> ¢(3/2 + k)(—9°/2)*/h = | 


k ) « @ 


¢(3/2 + s)I(—s)(2°/2)’ ds = 1, , (A2.4) 


the contour extending along the imaginary axis and including the pole at s = 0. Shifting 
the origin to the left [R(s) < 0] so as to include the single, simple pole of the ¢-function, 
here at s = —1/2, we easily show by Cauchy’s theorem that 


(8/2 + 9r(—s) & - (9°/2", 


i= (Residue at s = -}) . . 


J —~wi— 


(A2.5) 
(c > 1/2), 


since the contribution of the integrand along the infinite arcs (+ © i to + i — c) van- 
ishes.** Equation (A2.5) then allows us to write, 
*Whittaker and Watson, ibid., Ch. XIII. 
**For we note that 
| '(—s)(2?/2)*¢(3/2 + s) | = exp (—2 | v|/2)-| | 


as v->+©, where s = u + iv, independently of the value of c. 
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at least for the leading term as Q =, J, = (27)'’’/Q, a result also obtained by Rice 
[9] and Blachman [3]; (see also Appendix V of [1] for the other terms). 
The case of no limiting (rm —©) may be handled in the same way. The contour 
integral representation of the series (3.7) of [1] is now 
(1/2)’, 


bs Gath(om ¢ 17 =P (- 8 am + 2) 











1 f** cos xsl'(—s)I(s + 1/2)’ . ds 
ofl boa a he Eafe) exp [—97/(28 + 2)] & A2.6 
claw + Det po” hte +8 (A3.4) 
~ T'(—s)I'(s + 1/2) * ds 
= oe nl omm—mae Exp [— 0°/(2s + 2)] —, A2.7 
s Td/2 — ore + Die +)? oP L-F/As + Blan, —(A2.7) 
where we have used the relation sin rz = 2/I'(z)I'(1 — z) to transform the cosine term 
in (A2.6). The integral (A2.7) is convergent for all values of s = u + w whose real 
parts are finite, for by (A1.6) we have 
lim | P(=s)FGs + 1/2) exp [—9°/(2s + 2)] | | 1 
sri ri/2—slet+ et"? | Jo | 


The leading term is then found as in (A2.5) from the pole at s = —1/2, so that 





= 2)? exp [—7/(2 2 ae Ee 
FUME awa EM weep (—o). 
The rapid falling-off of the spectral intensity as 2 —©o when there is no limiting is 
attributable to the equally rapid decay [— exp (—*)] of the spectral intensity of the 
wave entering the discriminator. Roughly speaking, the number of beat-frequencies, 
and hence the spectral density, produced in the discriminator for very large 2 decreases 


lim 


2 m= 


in the same way as the original spectral intensity. 
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HEAT CONDUCTION IN A MELTING SOLID* 
BY 
H. G. LANDAU** 


Ballistic Research Laboratories, Aberdeen Proving Ground 


1. Introduction. There are a number of problems in heat conduction in which one 
material is transformed into another with generation or absorption of heat. Examples 
are the melting or freezing of a solid and the progress of a temperature-dependent 
chemical reaction through a solid. Some of these problems have been treated under 
various assumptions.’ The problem considered here is the melting of a solid when the 
liquid is removed immediately on formation. This problem seems to have been con- 
sidered before only by Soodak,’ who gave the steady state rate of melting and carried 
out numerical integration by finite differences for some special values. 

We first state a general problem and a few results. Then the problem is specialized 
to make possible a more detailed treatment. For the case of a semi-infinite melting® 
solid with its face heated at a constant rate, the exact steady state solution is derived, 
and the problem is put in a form involving only a single parameter. Numerical in- 
tegrations have been carried out for this case for several values of the parameter ex- 
tending over the entire possible range of the physical constants involved. 

2. Statement of problem. We consider one-dimensional heat conduction in a solid 
which initially extends from x = 0 to x = a. The face x = a is insulated, while heat 
flows in through the other face at a rate H(t) per unit area. If heating continues long 
enough, the face at x = 0 reaches the melting temperature 7',, and melting commences. 
The liquid is removed immediately on formation, say by being blown away, so that 
the face of the solid, which was initially at = 0, moves inward and at time ¢ is at the 
position x = s(t). If the heating rate decreases enough, the face temperature may drop 
below T’,, so melting stops temporarily and s(t) remains constant; then melting may 
resume if H(t) increases sufficiently. The temperature of the solid T(z, ¢) and the thick- 
ness melted s(¢) are the quantities we wish to find. 

In addition to the terms defined above we use the following: 

c = specific heat, p = density, k = thermal conductivity, x = distance from initial 
position of heated face, T(x) = initial temperature distribution, Z = latent heat of 
fusion. 

Then the equations describing the process will be as follows: 


ar _ 2 (x 22) . 
CP > = 5, k an)” s(t) <x <a, t > 0, (2.1) 


*Received Nov. 19, 1948. 

**Now with the Committee on Mathematical Biology, University of Chicago. 

Several references are given in H. S. Carslaw and J. C. Jaeger, Conduction of heat in solids, Oxford, 
1947, p. 71. An additional reference is L. I. Rubinstein, Isvestia Akad. Nauk SSSR, Sér. Géograph. 
Géophys. 11, 37-54 and 489-496 (1947). 

2H. Soodak, Effects of heat transfer between gases and solids, Ph.D. thesis, Duke Univ., Durham, 
N. C. (1943). 

*Although the term “melting” is used, the treatment in this paper would apply if the phase change 
were directly from solid to gas by sublimation. 
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T(x, 0) = T,(z) < T,. , e425 4 t= 0, (2.2) 
oT ; 
a 0, z=a, t> 0, (2.3) 
H(t) = —k _ + pL 4 2 = a(f), t> 0. (2.4) 


The last equation expresses the fact that the heat input H(t) equals the rate of heat 
flow into the solid plus the rate of heat absorption by melting. Condition (2.4) can be seen 
to hold both during melting and non-melting if we specify that at the heated face 
x = s(t) 


ds » @, for T(s(t), ) = T.. , (2.5) 
dt 

S a for T(s(t), t) < T 2.6 
dt - ’ Or s )s ) me ( ° )) 


The essential difficulty in the problem is in the determination of the unknown 
moving boundary s(t). The problem can be seen to be non-linear because two different 
heat input functions, H.(t) and H,(t), would give different boundaries, s,(¢) and s,(¢), 
so that the two solutions would not apply in the same region in the z,t plane, and the 
solution for the sum H,(t) + H,(t) cannot be obtained by addition. In this paper the 
existence and uniqueness of 7'(z, t) and s(t) are assumed. 

The moving boundary can be eliminated by the transformation § = (a — x)/[a — s(t)]. 
Then the boundaries are fixed at § = 0 and é = 1. Equations (2.1) to (2.6) become 


aT beat) 12 (2) 
Oreste =G—s a xe) 0O<§<1, t>0, 
T(x, 0) = T,(é) < T,. , 0<t< 1, t= 0, 
aT, 
="™ & = 0, i> G, 
k oT ds 
HO = ag te | 
| 
*> 0, for T(l, ) = Ta te = l, t>0 
dt | 
_ 0, fort, <T 
oC , »D m 


With the equations in this form the non-linearity of the problem can be seen from 
the fact that s and its derivative, which occur in the coefficients of the differential 
equation, depend on 07'/dé at = 1. We do not use the equations in the above form 
but they might be useful for numerical integration of the general problem. 
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Gauss’s theorem applied to the heat conduction equation (2.1) gives 


“9 | 2 (« or) ~ oF] dx dt = [ [x22 a+ cot ae, (2.7) 


where D is any region in the z,t-plane to which the theorem applies, and C is its boundary. 
Applying (2.7) to the region bounded by the lines t = 0, t = t, , x = a, and the curve 
zx = s(t), we obtain for its two terms: 


[era = [ | a — pb. | au = [He at ~ plate), 


[ ceT dx = , epT (x) dx — [ cpT (x, t;) dx — [ corto, t) de) dt 


“Cc /0 8(t,) /0 





In the last integral ds/dt = 0 except when T(s(é), t) = T,,, ; this integral, therefore, equals 
cpT ,,s(t,), and we have 


et 


[ H@at=p] (eT (e, t) — P(x) de 


(2.8) 
+ | tcf) — | ePo(z) az|, 


where we have dropped the subscript 1 from t, . This equation, as well as (2.7), holds 
when k and ¢ are functions of 7’. It can be seen that (2.8) equates the total heat inflow, 
up to time ¢, to the change in heat content of the remaining solid plus that of the liquid. 

Since T(x, t) < T,, , the first term on the right of (2.8) is bounded. Thus if {5 H(é) dt 
increases to large enough values, then s(t) must finally reach its largest possible value, 
s(t) = a; that is, if heating continues long enough the entire solid melts. The time, 


t, , when this occurs will be given by 


at 


[Hw at = aol L +P, — a? f eP(2) ae |, (2.9) 


“0 “0 


The rate of melting at this time can be obtained from (2.3), (2.4), and the con- 
tinuity of @7'/dx. The result is that at the time ¢, when melting is completed, 


ds(t,) _ H(t.) (2.10) 


dt pL i 


During melting, 97'/dx < 0 at x = s(t), otherwise T(z, t) > T,, for some x. Conse- 
quently from (2.4), ds/dt < H(t)/pL. Hence if H(é) is non-decreasing, the melting rate 
will reach a maximum at the time ¢, . 

The value of T(z, ¢) during any period when there is no melting is not difficult to 
find, under the assumption that c and k are constant. This value can be found by various 
well-known methods such as the use of the Laplace transformation or the method of 
sources and sinks.* The result is that before melting starts 


‘Carslaw and Jaeger, loc. cit., ch. X. 
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na vad - P 
T(x, t) = (4b’rt)"'” | T(t) Do {exp | -& Ae, — | 


Jo n=—@ 


. —— _\2 
—_ a si a (2.11) 





at ‘ , 
+ (akcp)~'”” | H(t — 7) b> exp | 25 pro | a 

Jo pentesl 4b°t T 
where b’ = k/cp is the thermal diffusivity. An equivalent expression in the form of a 
Fourier series could also be written. The time when melting commences, ¢ = ¢,, , will 
be given by T,, = 70, t,,). If melting stops at any later time ¢ = ¢, , then during the 
period when 7'(s(t,), t) < 7, , the temperature will be given by the expression similar 
to (2.11) but with the origin of distance and time coordinates at s(t,) and ¢, . 


Tg 2 y2 f° 1 Lo fe — t+ 2n(a — s(t,))]° 
T(x, t) = [4b°’(t — ¢,)]7'” | Té, t) > {ex E cr — & + én(a - L | 
(t, £) |4b°r )| T(é, t,) exp 451 i) 


“s(t,) n=—© 


Iz+t—2 2ns(t,)}? | 
+ exp B= on x or} dé (2.12) 
ey << [a — 2na + (2n — DaeoT | dr 
ke “ ap | —-2— SS ee fe". 
+ (rkcp) | H(t — 7) p>) exp | 45° —1) G— 4)” 


Instead of the heat input rate H(t) being given as a function of time, the problem 
could also be stated with heat input defined by the usual “radiation” condition as 
h([T,(t) — T(s(t), })], where T, is the temperature of the radiating medium. If T, is 
large compared with the surface temperature, the two statements will be practically 
equivalent. In the rest of this paper we take H as constant, which is the same as putting 
h and T, constant during melting. 

3. The semi-infinite solid. We now consider a much more special problem. In many 
cases occurring in practice the melting will proceed only a relatively small distance 
into the solid, so that it will be sufficiently accurate to take a = ~, that is, a semi- 
infinite solid. At the same time we assume that the thermal properties c and k, heat 
input rate H, and initial temperature 7’, are constant with H > 0. In some applications 
the variation of these quantities may be important. Even for these cases the results 
given here are useful as an approximation, and in showing how the progress of the 
melting depends on the physical quantities involved. By assuming these parameters 
to be constant we can also give a more detailed treatment which covers the entire range 
of physically possible values of the parameters. 

The condition (2.3) now becomes 


= x26 r—-@ i>0 (3.1) 


which is equivalent to 


T—T,, £-0,t>0. (3.2) 
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rhis equivalence can be seen from 


ar’ T ar 
I (2) Reding [ E de + T= at| (3.3) 
D 


~ 
4 


applied to a region D for which z is sufficiently large. Equation (3.3) like (2.7) is a form 
of Green’s theorem for the heat conduction equation,’ and can be obtained by applying 
Gauss’s theorem to 0(T”)/at — b’d*(T’)/dx” where T satisfies Eq. (2.1). 

The solution of the problem before melting starts is now, either from (2.11), or 


directly,” 


a t my 1/2 ( =) a (=:5) | 
T(x 9 j= Te + 2n(;) E exp Ab?t bi” erfc 2b”? , 
(3.4) 


t 1/2 2 
= T, 2 ) ierfe (=) 
a kep ed i 
where ierfe is the notation for the integral of the error function complement introduced 
by Hartree.’ 
From (3.4) the time when melting starts, ¢t,, , is given by 


— Thee ip _ = ( L ) 
tn = wiG a T.) =\p H bm} . (3.5) 
Here 
1/2 = 
ye En — Tr) 3.6) 


2 L 


is the ratio of the heat content change from 7, to T,, to the heat of fusion, multiplied 
by a convenient factor. By the transformations given below the whole problem during 
melting is reduced to a form in which the only parameter is m. 

In this semi-infinite case it is possible to eliminate the moving boundary by intro- 
ducing a new distance variable x — s(t), the distance being measured from the moving 
face of the melting solid.” Since we now have left to consider only the melting process, 
we shall measure time from t,, . We also make the variables dimensionless so as to obtain 
the equations in the desired simplified form. The new variables are 





. — Ty. ~ 

v= vz, y) = UT, — T.)’ (3.7) 
x— st 

z= = ae ) (3.8) 


5. Coursat, Cours d’analyse mathematique, vol. 3, Gauthier-Villars, Paris, 1915, p. 308. 
*R. V. Churchill, Modern operational mathematics in engineering, McGraw-Hill, New York, 1944, 


p. 107. 
7D. R. Hartree, Properties and applications of the repeated integrals of the error function, Mem. and 


Proc. Manchester Lit. and Phil. Soc. 80, 85-102 (1935-36). 
8Since x — s(t) = lime... a[l — (a — x)/(a — s(t))], this is a special case of the transformation defined 


by £ above, not an entirely different transformation for eliminating the moving boundary. 
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t 
aa ie l, (3.9) 
a pL 8 te 
o= i.” (3.10) 
_ do _ plds (3.11) 


athe dy A dt 


Then the equations for the melting process become 


Ov Ov , . Ow : 
. ay ay + myu(y) ay * @, y > 0, (3.12) 
v = ierte (2), a> 0, y = O, (3.13) 
v— 0 zZ—« y > 0, (3.14) 
Ov “ 
l= -2> + u(y), z= 0, y > 0, (3.15) 
Oz 
y= TF - Pie Z2= 0), y > 0. (3.16) 


This form for the equations is convenient because there is only one parameter m, 
and it occurs only in the differential equation (3.12); the initial and boundary conditions 
contain no parameters. This form for the equations, therefore, greatly simplifies the 
numerical integration which is described later. 

Mquation (3.15) should be regarded as the definition of the unknown reduced rate 
of melting, u(y). It is possible to eliminate u from the equations by substituting its value 
from (3.15) into (3.12). Equation (3.12) then becomes 


Ov a’v ‘ Ov(O, y) | Ov . 4 
= —; + m1 age 2] pee (3.12’) 


oy Oz Oz Oz 


which is clearly non-linear. 

4. Steady state solution. We first note that there is a steady state solution for the 
system (3.12) to (3.16); that is, the temperature of the melting, semi-infinite solid 
approaches a state characterized by the migration inward at a constant velocity of a 
fixed temperature distribution. It is found from the numerical integration that dv/dy > 0. 
Since v < x” it follows that dv/dy — 0 as y © and that v tends to a steady state. 
If we put dv/dy = 0, then yu is constant, say nu = u, , the constant steady state value 
of the reduced melting rate. From (3.12), (3.14), (3.15) and (3.16) this steady state is 
found to be given by 

v—v(z) = r '” exp (—mu,2), (4.1) 
with 


ue = (1 + 2mx'”)™. (4.2) 
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For the steady state rate of melting in terms of the original variables, this gives 





ds — H 
a’ "2st T.)|’ (4.3) 


which is also obvious physically. This does not give the steady state value for the thick- 
ness melted, but if we use, as a first approximation, a constant rate of melting, then 
the steady state thickness melted is approximately, 


o(y)=uy, or s(t) =V(t —¢,), (4.4) 
and the steady state temperature distribution in the original variables becomes 
T=T, + (T, — To) exp {—Vb-*[x — V(t — t,,)]}. (4.5) 


Actually the melting rate is initially less than its steady state value, so that (4.4) 
and (4.5) give values for the thickness melted and the temperature which are too large. 
We can obtain the exact steady state value as follows. 

Equation (2.8) written in terms of the reduced variables is 


o(y) = u(y +i- 2/ v(z, y) az). (4.6) 
Then by use of (4.1), the steady state value of o is 
a(y) > o.(y) = we(y + 1) — 2/n'*m, (4.7) 
or correspondingly 
s(t) — Vt — (k/H)(T,, — To), (4.8) 


and the exact steady state value for T is 
T—>T)+ (Tn — To) exp {—Vb"[xz — Vi + k(T,, — T.)/H)}. (4.9) 
Since v(z, y) increases with y, it can be seen from (4.6) that the steady state value 
gives a lower bound for the thickness melted. Also since 


[ v(z, 0) dz = [ ierfc (z/2) dz = 1/2, 


(4.4) gives an upper bound; that is, 


V(t — tn) > s(t) > Vt — (k/H)(T,, — T>). (4.10) 


The difference between these two bounds is 


KT, — To) _ # ( qs) | 
H [ 4 sa c(T,, a To F 


In the applications, this quantity is often small; then the steady state value is a good 


approximation of s(t) for all ¢ > ¢,, . 
By taking derivatives of (4.6) we find that the melting rate is always less than its 


steady state value: 


uly) << me, or & <¢ TF. (4.11) 
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io 2) 
1°) 


Physically, this is obvious since the heat content of the solid is increasing during the 
melting. 

5. The case m = 0. The parameter m in Eq. (3.12) can have any positive value. 
The value m = 0 does not correspond to a physically realizable situation but is a limit 
which is approached when the latent heat L becomes large. In this case, the system 
(3.12) to (3.16) can be solved in terms of tabulated functions. This solution is of value 
since it is one limit of the set of solutions for 0 < m <o. Also the solution for m = 0 
is an approximation to that for any finite m when y is sufficiently small, because u(0) = 0. 

For m = 0 the equations are linear. Equation (3.10) does not enter into the solution 
but is merely an equation for determining u(y) from the solution. The solution for 
m = 0 can be written as 


vo(z, y) = w'”” erfe (z/2y) 


+ (4ry)7'” [ ierfe (2) {exp | -¢ sf" — exp | -2+ 2 |} dg. (5.1) 


“0 





The indicated integrations are rather tedious; hence we omit the details and give only 
the final result. 
\ —-1/2 /é 2 —1/2) 
v(z,y) =7 erfe (2/2y) — (z/m) arc tan (y"'”*) 


+ [(y + 1)/x]'” exp [—2°/4(y + 1)] erf {z/2[y(y + 1)]'*} (5.2) 

+ 22 We/[2y + 1)]'’, y'”), 
where 

na abr , - 
Wa, 8) = (2m) | -[ exp [—(@* + y*)/2] dy de 
is a tabulated function.® 
This gives for the reduced rate of melting 
uo(y) = (2/m) arc tan (y'”), (5.3) 

and hence 

oy) = (2/m)[(y + 1) are tan (y'””) — y, (5.4) 


which as stated above are approximations to u(y) and o(y) for finite m when y is small. 
“a 


Since v,(z, y) can be used to start the numerical integration for 0 < m <o, it is 
useful to have simpler expressions than (5.2). For y > 0, vo(z, y) is analytic, and we 


can expand it in a power series in z. By use of duo/dy = 0°v)/dz’ and dv./dz 
1 J -1/2 . . . 
—7 are tan (y ) at z = O, this expansion is 
=) n 2n+1 
ee : a a ag 
v,(2,y) = 9? — —, (arc tan (y~"’")) =—— 
sn day ; (2n + 1)! 


(5.5) 


> 


_ _ 1/2 a, Seeere , —1/2 ee) —_ 
=f mS a tan (yj *) i2y'*y + + | y > @: 


°C. Nicholson, The probability integral for two variables, Biometrika 33, 59-72 (1943). Also, J. B- 
Rosser, Theory and application of {¢ exp (—z*) dx and fi exp (—p*y*) dy Ji exp (—z*) dz, Mapleton 
House, Brooklyn, 1948. 
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Also for z/y '” large 


Wey + 1IY2, y-) & F ef [e/2y + 1)" 


zx 1/2 
5, are tan (y’*), 


so that we get 
vo(z, y) = (y + 1)” ierfe [2/2(y + 1)”’]. (5.6) 


This relation shows that when melting begins, the temperature in the interior is nearly 
the same as it would be without melting. 

6. The case m = o. The case L = 0, which makes m = ©o, is of interest for its 
physical applications and also because it gives the other limit to the range of m. In 
fact for steel, with 7’) as room temperature, m = 27, which is “practically infinite” in 
the sense that the temperature and melting rates will be very close to those for m =~, 

For L = 0, the definitions of o and yu in (3.10) and (3.11) must be modified. We 
introduce the modified quantities 


2 Hs 
o= mo = ne? K(T,, = T.)’ (6.1) 
dw 7 cp ds 
= m = =” te (T,, — To) We (6.2) 
which are independent of L. In the system (3.12) to (3.16), Eq. (3.12) is now 
ov a°v Ov 
ay = ae +Myo, 2>0, y>Qd, (6.3) 


and (3.15) which came from (2.4) is now 


jn age 


a «2 = 9, oy > 0. (6.4) 


There is now no equation defining \, similar to (3.15) for u. However, from (3.16), 
dv/dy = 0 at z = 0, so that if the derivatives are continuous at z = 0, y > 0, Eq. (6.3) 


gives 
av Ov 
a2” + A(y) az = 0, .= 0, y > 0, 
or by use of (6.4) 
ap 
A(y) sal 2 ilies 0, y > 0, (6.5) 


which now takes the place of (3.15). The modified system of equations is then (6.3), 
(3.13), (3.14), (6.5), (3.16). 

For numerical integration this system is unsatisfactory because the determination 
of \ from (6.5) requires the calculation of a second derivative. Since this is not easily 
done with accuracy, we consider instead a derived system for dv/dz. 
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If we put 


u(z, y) = —2 5 4 ven? — ; [ u(t, y) df, (6.6) 
Oz Jo 


then u(z, y) will satisfy 


- = a + ry a z>0, y>QO, (6.7) 
u = erfc (z/2), > 6; y = 0, (6.8) 
u— 0, Zo, y = 0, (6.9) 

My) = —St, z= 0, y > 0, (6.10) 
u=tl1, z= 0, y > 0, (6.11) 


and it is easy to verify that if u satisfies (6.7) to (6.11) then v will satisfy (6.3), (3.13), 
(3.14), (6.5), (3.16). 

An important difference between the present case, L = 0, and L > 0, is that for 
L = 0 the melting rate is not zero at the start of melting; instead we now have 


(0) = 2” (6.12) 
corresponding to 
ds(t,,) 2H — 
ts ee (6.13) 
dt mcp(T', — T'o) 
For L > 0, ds (t,,)/dt = 0 was necessary for continuity of d7'/dz in t at ¢ = t,, , but with 
L = 0 this condition is not necessary, as can be seen from (2.4). 
The steady state solution is now instead of (4.1), 
v. = w  ” exp (—2'”*2/2), (6.14) 
and the steady state values of \ and w are 


A, = wr”? /2. and wo =rA.(y +1 — 4/n). (6.15) 


The steady state equations in terms of the original variables are unchanged except that 
now L = 0. 

7. Numerical integration. For m > 0, the solution was obtained by numerical inte- 
gration using finite differences. We describe the method used, but do not go into a dis- 
cussion of methods for numerical solution of the heat conduction equation and their 
convergence and stability. These questions have been discussed by several writers,” 
although much remains to be done, especially for non-linear problems. 


OR, Courant, K. Friedrichs and H. Lewy, Uber die Partiellen Differenzenglieichungen der Mathema- 
tischen Physik, Math. Ann. 100, 32-74 (1928); N. R. Eyres, D. R. Hartree, et al., The calculation of variable 
heat flow in solids, Proc. Roy. Soc. Lond. (A) 240, 1-57 (1946); C. M. Fowler, Analysis of numerical 
solutions of transient heat-flow problems, Q. Appl. Math. 3, 361-376 (1946); J. Crank and P. Nicolson, 
A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction 
type, Proc. Camb. Phil. Soc. 43, 50-67 (1947). Unpublished methods proposed by J. von Neumann. 
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We use Az and Ay for the intervals in z and y, and write v(i, j) for v(z, y) = v(iAz, jAy) 
and yu; for u(jAy). The time derivative is approximated by a forward difference ratio 
and the space derivatives by central difference ratios, 


2b D ew (ay) fli, f+ 1) — v66, (7.1) 
OND mw (202) [oli + 1, f) — 0 — 1H) (7.2) 
OG D ew (a2)-* [uli + 1, f) — QE, ) +06 — 1, H], (7.3) 


Oz” 


so that Eq. (3.12) is approximated by the difference equation 


— nail At 
v(t, 7 + 1) — oft, 3) = - ¥ 


=i {ue + 1, j) — 2v(i, j) 


(7.4) 
+ v(t — 1,9) + mp; 2 [ot + 1, 7) — ot — 1 a}. 


It is known that the solution of the difference approximation to the heat conduction 
equation will not be stable unless the ratio Ay/Az’ < 1/2. Stability means that errors 
which occur in the process of solution, such as rounding-off errors, are damped out 
instead of being magnified. We then take Ay/Az’ = 1/2, which simplifies (7.4) by re- 
moving the term in v(7, 7), so that the difference equation becomes 


v(i, j + 1) = $1 + mp,Az/2)06 + 1, J) 
(7.5) 


+ (1 — my, Az/2)v(¢ — 1, J). 


To compute yu, a third order formula for dv/dz at z = 0 was used. Then from (3.15) 
and (3.16), 


uy = 1 + (Az)""[—11/32'” + Gv(1, 7) — 3v(2, 7) + 30(3, J)]. (7.6) 


The computation proceeded along “‘lines,’”’ the values for fixed j and all 7 being computed 
before going on to j + 1. The initial values are given by (3.13), or we can start with 
j = 1 using the solution for m = 0, Eq. (5.2). Having v(?, 7), (7.6) is used to find yu; 
and then (7.5) is used to compute v(z7, 7 + 1). The computation was carried to values 
of i large enough to make v(i, j) < 5 X 10°’. 

The calculations were carried out partly on the Bell Relay Computer and mainly 
on the I B M Relay Calculator"’ at the Ballistic Research Laboratories for m = 0.2, 
1, 2,5, 10 and . The modifications necessary for m = © are clear from Sec. 6. 

In the applications of the results a very high accuracy in uw and a is not needed. 
The calculations were carried out so as to give uw (which is the quantity with the largest 


''These machines are described in F. L. Alt, A Bell Telephone Laboratories computing machine, Math. 
Tables and Aids Computation 3, 1-13, 69-84 (1948), and W. J. Eckert, The I B M pluggable sequence 
relay calculator, Math. Tables and Aids Computation 3, 149-161 (1948). 
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error) generally correct to three significant figures, although the third figure was some- 
times in doubt. The accuracy was checked in various ways, especially by using Eq. (4.6). 

Since » and v change rapidly at first, it was necessary to start the computation with 
small values of Az and Ay. When the rate of change decreased, the interval size had to be 
increased to keep the computation from becoming too lengthy. As the steady state 
solution was approached, it would have been desirable to use relatively large intervals 
of reduced time Ay, but the limitation Ay/Az* < 1/2 would then have required Az to 
be so large as to introduce serious errors. The limitation on the relative size of the time 
and the space interval could have been avoided by using a more complicated difference 
approximation. 

To extend the values of » to the steady state value, a semi-empirical method was 
used, as follows. If in Eq. (3.12), mu(y) is replaced by a constant r, then the system (3.12), 
(3.13), (3.14) and (3.16) can be solved by using the substitution v = w(z, y) exp [—rz/2 — 
r’y/4], to eliminate the term r dv/dz. The solution is then, 


*X —rz/2 — ” } eae . /2 z2-—-¢ 4 
v(z, y) = “ae Ta y)? oo. | lerfe (¢ 2)e"*” {exp | _{ on | 
= €Ap | - o” I dt ( 
. by J . 
+ = 


l" io: gE 2? | dn 
“ e) —_ ie 
2r Jo 4 4(y — nj (y — 0)” 


The second term v,(z, y) can be evaluated as 


ve(z, y) = (4x)7'”” | ext (: aa) +e” erfe (=) | (7.8) 
2y 2y 


“J 
~J 


As y >, the first term in (7.7) approaches zero, while v,(z, y) ~ 7 ““e "*, the steady 
state solution. We then use dv.(0, y)/dz to indicate the dependence of u on y, for y large 


enough so that u is almost equal to u, . From (7.8) 


ov.(0, /) ‘ 7, 1/2: . /2 my 
er ~ Ys —x" tir + y'” ierfe (ry'’”/2)], (7.9) 
Oz 
Putting r = mu, and using (3.15) we have 


w= u.[1 — 2(ry)”'” ierfe (mpy'’”*/2)], 


(7.10) 
> u[1 — 2(ry)'” ierfe (mp.y'”*/2)], 


for large y. 

It was found that for y large the computed values of uw varied with y in this way, 
except that the second term in (7.10), which is a small correction term, had to be multi- 
plied by a constant depending on m. Thus yu could be extended to arbitrarily large y 
with an error which is within the desired limits, because the extension was started for 


uw within a few percent of y, . 
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Figures 1 and 2 show the values of the melting rate plotted against time for m = 0, 
0.2, 1, 2, 5, 10 and . In Fig. 1 it can be seen that for all m, the values of « approximate 


10 


6 
a * | 
a|x * 
u 





01 02 03.04 .06 08! A a4 828% . 2 
yee -I 
m 
Fig. 1. Rate of melting of the semi-infinite solid for various values of m: 
m = r'2¢(T'm — T'o)/2L, ds/dt = rate of melting, p = density, L = 
latent heat, H = heat input rate per unit area, = time, tm = time 
when melting begins. 
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Fic. 2. Rate of melting as fraction of the steady state rate: 
V = H/p([L + c(T'm — To)], steady state rate of melting. 


those for m = 0 for small values of y, and that u approaches its steady state value more 
rapidly the larger the value of m. In Fig. 2 the ordinate is u/u, (for m = @, it is \/A,) 
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so that all the curves approach unity as an asymptote. In this figure it can be seen that 
as m becomes large u/u, near y = 0 approaches the value for m = ~, although the value 
of u/u, at y = 0 is quite different from that of m =. Figure 3 shows the thickness 
melted; the asymptotic values here are ¢/u, > y + 1 — 2/(x'”muy,) for m > 0, and 
y+1— 4y'’/x for m = 0. 





-O1 02 .03 04 .06 .08 |I a 3 4 6 8 10 2 


Fic. 3. Thickness melted, s. 
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—NOTES— 


A NOTE ON RIEMANN’S METHOD APPLIED TO THE DIFFUSION EQUATION* 
By JOHN W. MILES (University of California at Los Angeles) 


Summary. The method of Riemann (or Green) is applied to the integration of the 
parabolic equation in m dimensions and two variables. The result is applied to the two 
phase diffusion problem. The reduction to tractable equations is not simple, however, 
and the results of the investigation appear to be essentially negative. 

1. Introduction. The method of Riemann for integrating partial differential equa- 
tions has received considerable attention in the case of the hyperbolic equation”’®, 
where it has proved to be a most elegant approach. Its application to equations of 
the parabolic type has, on the other hand, received considerably less attention.** 

Pascal’ has considered the special case of fixed boundaries (in space), while Rade- 
macher and Rothe*® have considered the more general one-dimensional case where the 
dependent variable is specified over an arbitrary boundary but did not apply it to the 
solution of particular problems. The following note will be restricted to the case of 
two independent variables, but not necessarily to one dimension. The extension to an 
arbitrary number of independent variables may be made along the lines indicated by 
Webster.” 

2. Riemann integration. The equation to be integrated is 


O | m1 O gt 9 ie a 
L(u) = = E aa u(x, 0 | Kr ay Ue t) Q(x, 2). (1) 


If (1) is regarded as the equation of heat conduction, then u(z, ¢) is the temperature, 
x is a space coordinate, ¢ is the time coordinate, x is the diffusion constant, and Q is 
the source strength. The space is m dimensional, and corresponds to one dimension, 


*Received May 17, 1948. 
1H. Bateman, Partial differential equations of mathematical physics, Dover Press, New York, 1943, 
p. 126 ff. 

2R. Courant and D. Hilbert, Methoden der mathematischen Physik, J. Springer, Berlin, 1936, vol. 2, 
p. 311 ff. 

8A. G. Webster, Partial differential equations of mathematical physics, G. E. Stechert and Co., New 
York, 1927, p. 244 ff. 

‘P, Frank and R. v. Mises, Differentialgleichungen der Physik, Friedr. Vieweg & Sohn, Braunschweig, 
1925 (7th Ed.), vol. 1, p. 605 ff. 
5E. Pascal, Repertorium der Héheren Analysis, B. G. Teubner, Leipzig, 1929, vol. 3, pp. 1175-1187. 
6. J. B. Goursat, Cours d’analyse, Gauthier-Villars, Paris, 1935, vol. 3, p. 308 ff. 
**References (5)-(8) were brought to the author’s attention by the reviewer, who also pointed out 


R 
that the method under discussion is designated as Green’s method therein. It is the author’s impression 
that Green applied the method only to the elliptic equation, and that Riemann developed the generaliza- 
tion to other equations and to an arbitrary number of variables, cf. the discussion given by Webster, 
loc. cit. (3), p. 239. See also J. Hadamard, Lectures on Cauchy’s problem in linear partial differential 
equations, Yale Univ. Press, New Haven, 1923, p. 57. 

7Pascal, loc. cit., p. 1178. 

8Frank and v. Mises, loc. cit., chapter by H. Rademacher and E. Rothe, p. 646 ff. 

*Webster, loc. cit., p. 257 ff. 
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two dimensions with axial symmetry, or three dimensions with spherical symmetry for 


m = 1, 2, or 3, respectively. 
The operator which is adjoint to L is defined by 
0 = 0 a 
M(v) = - x” '—vu(z, 0| + Kx ' — o(z, t). (2) 
Ox Ox ot 
If the combination vL(u) — uM(v) is integrated over a surface S in the (é, r) plane, 


which is bounded by the curve C, it can be shown that” 


' wv au he 
// (vL(u) — uM(v)] dé dr = p { u = —v ] cos (n, &) + xuv cos (n, be ‘dl = (3) 
: og dg J 
8 


A 
where (£, 7) are running (2, ¢) coordinates, and cos (n, £) and cos (n, 7) are the direction 
cosines of the inwardly directed normal to C with respect to the & and 7 axes. 


, 
wt) “© 





g 








Fia. 1. Curve in (¢, r) plane for integration of Eq. (3). 


The identity (3) will now be applied to a curve consisting of two segments C, and 
C, , as shown in Fig. 1. C, is a portion of the line r = ¢, and C, is a curve for which 


7 < t, but is otherwise arbitrary. Moreover, it will be assumed that 


M(v(é, r)] = 0, (4) 
xt” ‘v(é, t) = (x — &), (5) 
where 
d(x — §) = 0, xrHFE (6a) 
| d(a — §) dt = 1, e>0 (6b) 





This follows directly from the more general result 





derived by Webster, loc. cit., p. 245. 
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i.e., v(~, 7) is a solution to the adjoint equation (4) which reduces (essentially) to the 
Dirac delta function, as defined by (6), at 7 = t. Now along the curve C, 


dl cos (n, ~) = 0, (7a) 


dl cos (n, r) = —dé, (7b) 


‘ 4» 4 
? {[ > —v ou | cos (n, —) + xuv cos (n, oe dl 


=— p ul, t) d(x — &) dé, (7¢) 
= —u(z, t). 
Substituting L(w) = —Q from (1), M(v) = O from (4), and the line integral along C, 
from (7c), (3) reduces to 
ulz, =) = i/ Q(E, tvo(a, t; —, 7) dé dr 
(8) 


+ p {[: > — u 2 dr + xuv age 


where the adjoint function v is now regarded as a function of both the parameters 
(x, ¢) and the running coordinates (&, 7). 

3. Comparison with the method of characteristics. The foregoing solution is a 
degenerate case of the method of characteristics, as applied, e.g., by Volterra to the 
wave equation.'’ The two characteristic lines utilized in the solution of the hyperbolic 
equation in two coordinates merge into a single characteristic for the parabolic equa- 
tion. Moreover, in applying the method of characteristics to the hyperbolic equation, 
the curve C, is restricted to intersecting a given characteristic only once, whereas in 
the present case the curve necessarily intersects the characteristic twice. It may also 
be remarked fhat the adjoint functions appropriate to the hyperbolic equation are 
generally singular along the entire characteristics through the point at which the solu- 
tion is desired, whereas in the parabolic equation the adjoint function is singular only 
at that point. 

The latter property is perhaps more closely related to the solution of the elliptic 
equation, where the adjoint function is singular at a point, although this point is then 
located within C, rather than on its boundary. 

4. Boundary conditions. It will be assumed that the boundary condition is of the 
open Dirichlet-Neumann (mixed) type, viz. 


u(é, tr) + B(E, 7) < u(t, r) = &” "f(é, 7), (—, r) on C, . (9) 


\Webster, loc. cit., p. 266 ff. 
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If 8 = 0 the boundary condition is simply of the open Dirichlet type, whereas if 8B = © 
(with f replaced by B°'f), it is of the open Neumann type. (The adjective “open” 
implies that boundary conditions are not applied over a closed curve in the &, 7 plane, 
i.e., the curve is open in the direction of positive time.) Substituting the condition (9) 
in Eq. (8) to eliminate (du/d£), it is found that the explicit presence of u in the inte- 
grand may be eliminated by requiring the adjoint function to satisfy the boundary 


condition. 


a 
u(x, t; &, 7) cos (n, 7) + E v(x, tt, r) + B '(E, r)v(£, | cos (n, —) = 0, 


¥: 
= OS (10) 
(tf, 7) on C, 
whence Eq. (8) reduces to the alternative forms 
u(x, t) = [| QE, r)v(a, t; £, 1) dé dr 4 $ 8 (gE, r)f(E, r)v(a, t; §, 7) dr, (11a) 
[| QE, r)v(a, t; &, 7) dé dr 
(11b) 
af y i. 
c 
> f(E, r)| wv(x, t; &, 7) dE — — (2, t; , 7) dr |. 
$16 o[ ate te, nae Ze ud 00] 


Cs 


The result (11) replaces the problem of solving Eq. (1), subject to the inhomogeneous 
boundary condition (9), with the problem of solving the adjoint equation (4), subject 
to the homogeneous boundary condition (10) and the singular condition (5). It re- 
mains to show that the latter solution exists and is unique. The proof is evidently the 
same as that for Eq. (1), if the sign of ¢ is changed, and it follows that a necessary 
condition is that the boundary condition on v be open in the negative r direction (note 
that (5) closes the boundary condition in the positive 7 direction). Now, since cos (n, r) 


and cos (n, £) cannot vanish simultaneously, the only conditions for which Eq. (10) is 


indeterminate are 


cos (n, ) = 0] 


(t, r) on C,. (12) 


g, 


Bg, 7) = 0) 


Hence, for the smallest value of r (the initial time) on C, , assuming the curve to be 
continuous so that cos (n, ¢) = 0 at the minimum, 8 must vanish, and it is therefore 
necessary, ef. Eq. (9), that u be explicitly stated there. The physical implication is that 
the phenomenon described must have an initially prescribed state. 

5. Reduction for fixed boundaries. Perhaps the simplest contour C, consists of two 
lines of constant ¢ plus a line of constant r(<¢). For convenience, the last line may be 
0 with no loss of generality. The first two lines will be chosen as § = 2, 


taken as r 
0, would be a rather special case 


and ¢ Xx» , aS shown in Fig. 2. (Choosing the line & 


for m > 1.) 
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The boundary conditions (9) and (12) reduce to 


td, ey, 098 : u(z,,) =ai"f(), t= 1,2 (13) 
r = 0: u(t, 0) = w(t) (14) 

while Eq. (10) reduces to 
Wt, bt, +B. Sve, Ga, 7) =0, i= 1,2 (15) 

and Eq. (11b) becomes 
uz, ) = [ dr [  QG, Aula, 58, 0) ae 
+ xf wolGo(a, #5, O8"* ae 

(16) 


-}- [ Si(r) eM t;x2,, 7) dr 


“ [ f(r) ra t;x., 7) dr 


(xt) 














é 





E=x, U=U, E=X> 


Fic. 2. Curve in (£, r) plane for fixed boundaries. 


The result (16) is generally deduced directly from Green’s theorem, rather than the 
more general identity of Eq. (3). 

6. Application to two phase problems. As an example of a more general type of 
problem to which the result (11) could be applied, the freezing or melting of a solid 
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will be considered. It will be assumed that the substance undergoing such a change of 
state is located between two fixed boundaries § = x, , x. . Between these two boundaries 
lies a moving boundary, designated 2,.(r), as shown in Fig. 3. The temperature on 
the two sides of this boundary will be designated as u, and u, . 


. 





m 
" 


U, U= Us 
&=X,2(T) 








z 








=X, U=U, E=X2 
Fic. 3. Curve in (¢, 7) plane for two phase problems. 
At the boundaries § = xz, , § = zx, , and tr = O the boundary conditions (13) and 
(14) are appropriate, while the boundary conditions at z,, may be written 


(Zig , €) = Uee(Xig , 2) = 0, (17) 


ee oe 2 | datie(t) 
E a U(%2,t) — k re U2(Xi2 , 0 | =p ~ (18) 


where (17) states that the boundary is one at which freezing occurs, and (18) follows 
from a heat balance of a small element including a section of the boundary, k being 


the conductivity, p the density, and \ the latent heat. The result (8) may be applied 
separately to the regions to the right and left of the boundary z,, , where 


dl,, = ¥dr sec (n, &,,). (19) 
The adjoint functions in the two regions satisfy Eq. (5) at r = t, Eq. (14) at & = a, 
and ¢ = z,, and 
v,(z, t;%i2,7) = v(x, t} 2, 7) = O (20) 
Substituting the conditions (13)-(15), (17), and (20) in Eq. (11b) and assuming 
Q = 0, the temperatures are found to be 


eZ12(0) 7 at ; F,) 
“u(¢, 0) = &k | U(E)v, (a, t; E, OE” dé + | Fi(7) ag v,(z, t; x, , 7) dr, (21a) 
nts at a 
U(x, t) = Ko | Uo(E)v2(x, é, oe" dé — fl T) aE v(x, * Le ) T) dr. (21b) 
* 0) “0 Ss 


Zia 
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The determining equation for x,. is now obtained by substituting Eqs. (21) in Eq. (18) 
with the result 


— 0 Za 
hses J ult ” 158 OR dE — Kane | o) & vvlare » t5 &, OR" ak 
ie 212(0) 
k, [ fil _% Ui(X2 , b; ny, a) drt he [fal fol Dae atlas t X., 7) dr (22) 


d 
= Xp dt X2(t). 


The solution of Eq. (22) for x,,(¢) will evidently be quite complex in practical appli- 
cations; moreover, the function v, and v, must be determined first.* For m = 1, 2, =@, 
and u(2,;) = const., a solution can be given and was first obtained by Riemann,” using 
similarity solutions to the original partial differential equation (1). The physical problem 
investigated by Riemann was the freezing of a deep lake for constant surface tempera- 
ture, and his approach yields the desired results with considerably less trouble than 
the solution of Eqs. (21), (22). 

The problem of the growth of ice on cylindrical cables or spherical shells and the 
freezing of flat, cylindrical, or spherical ingots are all problems of great practical im- 
portance which fall in the category under consideration. Only the first of these problems 
has (to the author’s knowledge) been treated analytically, and then only quite ap- 
proximately.'* The author has attempted to solve such problems through the formula- 
tion of Eqs. (21) and (22) but was unable to obtain anything approaching an analytical 
solution in closed form, except for the simple problem of freezing a slab with a (required) 
constant temperature gradient at the interface of freezing. 


Perhaps the most expedient approximate approach is to assume an equation, say 
(1) 


(0 


x}> (t), for the propagation of the interface, introduce it in Eqs. (21) and solve for u; 
and u;'’, substitute these in Eq. (18) to obtain the next approximation x\2 (t), solve 
for u,”’ and u;”’, ete. While the convergence of this process would appear to be quite 
rapid, the evaluation of the integrals in Eqs. (21) remains difficult (even when v, and 
v. are readily determinable), and it will generally be necessary to have recourse to 
numerical or graphical techniques. 

7. Conclusion. While the method of Riemann brings about a rather general solution 
to the parabolic equation, it does not appear to present an immediately practical method 
of solving problems which are more complex than those which are readily treated by 
more standard methods. The general solution presented herein may nevertheless prove 
useful in leading to various approximate solutions, or it may be that the equations such 
as those presented herein may actually be solved analytically or on automatic com- 


puting machines. 





*In general, they will be analogous to the Green’s functions for fixed boundaries. 
2Frank and v. Mises, loc. cit, vol. 2, chapter by Riemann, 220-222. 
13C, L. Perkins and L. B. Slichter, Problem of ice formation, J. Appl. Phys. 10, 135 (1939). 
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ON FREE VIBRATIONS WITH AMPLITUDINAL LIMITS* 
By AUREL WINTNER (The Johns Hopkins University) 
1. Consider a differential equation of the form 
y + a,(t)y + a(t)y = 0 (1) 
or, more generally (placing 
m=y May (2) 


and n = 2), a system of n differential equations 
‘ n 
t= > aul(t)rr , (@ = 1,---,n), (3) 
k=1 


where every coefficient function a(t) is given as continuous for large positive t, say for 
i’ < t <o. If A denotes the matrix (a;,), and z the vector (x, , --- , z,), then (3) can 
simply be written as 


z = A(t)z. (4) 


Since n can be replaced by 2n, there is no loss of generality in assuming that a,, , xz; in 
A, x are real-valued. This will always be assumed in what follows. 

The trivial solution, x(t) = 0, of (4) will be excluded. Then, if x(¢) is any solution 
vector of (4), there cannot exist any ¢, for which x(t.) = 0 (where 0 is the zero vector). 
For, on the one hand, z(t) = 0 is a solution of (4) satisfying the initial condition z(t.) = 0 
and, on the other hand, any initial condition determines a solution x(t) of (4) uniquely. 
Accordingly, if r(¢) denotes the length of the vector x(t), then 


a(t) = r(t)e(t), where r(t) > O and | e(t) | = 1 (r = |2z/)). (5) 


The positive scalar r(t) and the unit vector e(t) will respectively be referred to as the 
amplitude and phase factor of the solution 2x(é). 

Various conditions are known which, when satisfied by the matrix function A(é), 
will assure for the system (4) the following type of asymptotic behaviort: Every non- 
trivial solution vector z(t) of (4) tends, as t ~©, to a non-vanishing limit vector; in 
other words, every amplitude r(t) tends to a finite, positive limit r(), and every phase 
factor e(t) to a certain unit vector e(). From the point of view of the theory of vib- 
rations, there are two objections to any criterion of this type. 

The first objection is that, although such a criterion can be made applicable by 
using the method of the variation of constants**, it cannot apply directly to simplest 
vibration problems. In fact, even if the problem is that given by the case a,(t) = 0, 
a,(t) = 1 of (1), that is, by the linear oscillator y + y = 0, it is seen from (2), where 
x, = ccos (t — y), z, = —csin (¢ — y) in the present case, that the above phase factor 
e(t) fails to tend to a limit; although the limit, r(@), of the amplitude exists and is 
positive, since r(t) = | c| = const. > O (unless x(t) = 0). 


*Received Sept. 20, 1949. 
TA. Wintner, On linear asymptotic equilibria, Amer. J. Math. 71, 853-858 (1949). 
**A. Wintner, Small perturbations, Amer. J. Math. 67, 417-430 (1945). 
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The second objection is that any criterion of the type in question becomes too severe 
by necessity. For, if only r = | x], rather than x = re (that is, r and e) is desired, then 
there becomes involved that issue which, on the one hand, is precisely the hard part 
of the problem and which, on the other hand, was not required at all. In fact, if the phase 
factor of a solution is known, then the amplitude of the solution (hence, the solution 
itself) can be obtained by a quadrature. f 

In order to see this, it is sufficient to observe that, since r? = 2.2, hence rT = 2.2, 
scalar multiplication of (4) by x gives rr = x.A(t)x. It follows therefore from (5) that 
rr = r’e.A(t)e, hence (log r)’ = e.A(t)e. Consequently, if e = e(t) is known, then log 
r(t) follows by a quadrature. 

2. In what follows, a criterion will be developed which deals only with the problem 
of asymptotic amplitudes, without any reference to phase factors, and is therefore free 
of the above objections. In other words, a condition will be deduced which is not violated 
by vibration problems (such as y + y = 0) and which, when satisfied by the coefficient 
matrix A(t), is sufficient to ensure for the corresponding system (4) the following prop- 
erty: 

The amplitude function, r(é) = | x(t) |, of every non-trivial solution vector, x = x(t), 
of (4) tends to a finite, non-vanishing limit, as t ~. For the sake of brevity, such a 
system (4) will be said to be of type (*). 

If a prime denotes the operation of transposing a matrix A = (a,,), that is, if A’ = 
(a,,), then the criterion in question can be formulated as follows: 

Let \ denote the least, and yu the greatest, characteristic number (eigenvalue) of the sym- 
metric matrix (1/2)(A + A’), where A is any real matrix. For large positive t, say for r< 
t <<, let A = A(t) be a continuous function, and suppose that the corresponding continuous 
functions \ = X(t), u(t) are integrable over the half-line t° < t < @, that is, that the integrals 


T T 
lim A(t) dt, lim u(t) dt (6) 
T+@ T+@ 
are convergent. Then the system (4) is of type (*). 
It is worth emphasizing that f° | \(é) | dt = © or J” | u(t) dt = © is allowed, that is, 
that the absolute convergence of the integrals (6) is not required. 


3. First, if = (& , --- , &,) is any vector, then 
f.AE = 46.(A + ADE (7) 
is an identity, since 
EAE = DOD aukée and 38.(A + ADE = DOD Hau + avde& , 


where A = (a,,), A’ = (a;). On the other hand, if \ denotes the least, and u the greatest, 
eigenvalue of a real, symmetric matrix B, then \ is the minimum, and » the maximum, 


attained by the (real) quadratic form ¢.Bé on the unit sphere, || = 1. If this fact is 
applied to the matrix B = (1/2)(A + A’), it follows from (7) that 
A < &AE S wif | E| = 1. (8) 


Next, if r = r(t) denotes the amplitude, and e = e(t) the phase factor, of an arbitrary 
non-trivial solution vector x = x(t) of (4), then, as verified at the end of Section 1, the 
logarithmic derivative of r(t) is identical with e(t).Ae(t), where A = A(t). Since the 
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hypothesis, | ¢ |, of the inequalities (8) is satisfied by the vector ¢ e(t), it now 


follows from (8) that 


A(t) < d log r(t)/dt p(t). (9) 
Finally, if ¢° u <v, integration of (9) between ft uw and ¢ v gives 
/ A(t) dt < log r(v) log r(u) < [aco at. 
On the other hand, the convergence of the improper integrals (6) means that 
/ A(t) dt — 0 and [ u)dt-O0fuowv,v >, 
jut the last two formula lines imply that log r(v) log r(u) ~ 0 as uo, vom, 


This means that the logarithm of r() tends to a finite limit as -—> ©. Since this is equva 


lent to the statement that r(d) itself tends to a finite non-vanishing limit, the proof is 


complete 


A GENERALIZATION OF ALFREY’S THEOREM FOR VISCO-ELASTIC MEDIA* 


By H.S. TSIEN (California Institute of Technology) 


1. Introduction. lor the non-homogeneous stresses in isotropic incompressible visco- 
elastic media characterized by linear relations between the components of stress, strain 
and their derivatives with respect to time, T. Alfrey has shown (Ref. 1) that in the 
boundary value problem, the stress distribution is identical with that 


A 


case of the first 
in an incompressible elastic material under the same instantaneous surface forces 
similar result was obtained for the second boundary value problem where the displace 
ments at the boundary are specified. It is the purpose of the present note to generalize 
this theorem to isotropic compressible media for problems involving body forces. Only 
the first boundary value problem will be discussed, as the corresponding theorem on the 


second boundary value problem is self-evident 


2. First boundary value problem. Let the displacements along the 2, y, z directions 
* of the six strain components can be written as: 


be u, v, w. Then the typical expressions* 
Ou 
€, - , 
OX 
(1) 
Ou Ov 
1 — 
OY Or 
If the six stress components are denoted by o, , 0, , 6, , Try » Tys » Tez » the components 


*Received Sept. 7, 1949 
**Throughout this note, only typical expressions are explicitly given; other expressions can be readily 


obtained by cyclic permutations. 
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of body force by X, Y, Z, and the surface force per unit area by i Y, Z, the equations 


of equilibrium are 


Oo, OTs, OF ss , 

Ox Oy + Oz a= (2) 
Hlere the body forces X, Y, Z are the result of external field or agent and will not be 
identified with the inertia forces of the material. The inertia forces are here considered 
to be negligible as is actually the case for a wide class of problems. If 1, m, n are the 
direction cosines of the normal to the surface, then the surface conditions are: 


X = lo, + mr,, + n7,; (3) 


To determine the stress distribution completely, there are in addition six equations 
of compatibility: 
O'€, d’e, _ d’y,, eo 0’, 0 ( OV ye OV ex 1) (4) 
oy Ox” Ox dy’ ” dy dz Ox Ox Oy Oz 
It remains to specify the relations between the components of stress, strain and 
their time derivatives. These relations will be assumed to be linear, corresponding to 
problems of small strain. If in addition the material is assumed to be isotropic, then 
purely on the ground of invariance under space coordinate transformation it can be 
shown that the required relations have to be of the following form: 


Po, = Q(r\e + 2ue,) 
(5) 


Pte, = Quysy ’ 


where » and X\ are constants, and 


c es + €y + €, (6) 
The operators P and Q are time operators defined as: 
9” am—t 
P= —+4a,.,- ~ forts + dy 
at” ot” 
(7) 
a" a” 
O° 5 +4. 554" +h 
dl at" ' 


The a’s and b’s define the characteristics of the material. They could be functions of 
time, but not functions of the space variables. Thus a material with changing properties, 
caused by the drift’ towards thermodynamic and chemical equilibrium, can be also 
represented by these operators. 
By eliminating the strains between the compatibility equations (4) and the stress- 
strain relations (5), one has: 
2 2\ + Qu a°O d (= ay a5) 42 ox 
i. 3X + Qu dx” A + 2p \dx Oy Oz Ox 


™ 2+ 22 30 fs ar) 
¥ = 
P| Tov + 3A + Qu dx dy + Oy + Ox 


0 
(8) 
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where 
0@=¢0,.+ 90, 1+ 2; (9) 
The equations (8) are sufficient to solve the first boundary value problem: The surface 
forces X, Y, Z and body forces X, Y, Z are specified for all values of ¢. For any given ¢, 
these forces must be in equilibrium. The problem is to determine the distribution of 
stresses fulfilling these boundary conditions. 
Now let 

X = X*9g(2), Y = Y*9(0), Z = Z* g(t) (10) 
The starred quantities are functions of space coordinates only. Then for equilibrium the 
body force must vary with time in a similar way. Thus: 


X = X*g(t), Y=Y*%*g(t), Z= Z*g(t) (11) 


The stress components can now be written in the same manner: 


= (12) 


a, = o*9(t), Try = T2,9(t), 


By substituting equations (11) and (12) into equations (8), it is easily shown that the 
starred quantities satisfy the stress equations for a purely elastic medium with Lamé’s 
constants \ and u. By substituting equations (10) and (12) into the boundary conditions 
(3), it is seen that the starred quantities also satisfy their corresponding boundary 
conditions. Therefore, in the case of the first boundary value problem, the stress dis- 
tribution is identical with that in purely elastic material under the same instantaneous 
surface forces and body forces. 

To determine the displacements u, v and w, one introduces the unknown time function 
h(t) such that: 


u= u*h(t), v = v*h(t), w= w*h(t) (13) 


where the starred quantities are again functions of space variables only. When equations 
(13) are substituted into equations (5), u*, v*, w* are found to be the displacements of 
a purely elastic medium under the loading X*, Y*, Z* and X*, Y*, Z*. Furthermore, 


h(t) is determined by: 
Qh(t) = Pg(dt) (14) 


with the initial condition that at ¢ = 0, h and its first (n — 1) derivatives vanish. The 
function A(t) is thus universal in the sense that it depends only on g(t) and the properties 
of the material. The other characteristics of the problem do not enter into its determina- 
tion. In particular, the function h(t) may be directly determined experimentally on a 
pure tension bar with the tension varied with time according to g(t). 

By superimposing solutions, the time dependence of the applied forces can be gen- 


eralized as shown by Alfrey (Ref. 1). 





IT. Alfrey: Non-homogeneous stresses in visco-elastic media, Q. Appl. Math. 2, 113-119 (1944). 
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NOTE ON THE KINEMATICS OF PLANE VISCOUS MOTION* 
By J. L. SYNGE (Dublin Institute for Advanced Studies) 


In 1911 G. Hamel (Géttinger Nachrichten, Math.-Phys. Kl. 1911, 261-270) obtained 
an interesting result, which may be stated as follows: 

Let # be a finite plane region with boundary B. Then the equation AY = F possesses a 
solution y for which both y and dy/dn vanish on B if and only if F satisfies 


[ FU dx dy = 0, (1) 


U being an arbitrary harmonic function. 
In other words, for the existence of a solution with this double boundary condition, 
it is necessary and sufficient that the function F be orthogonal to the linear space of har- 


monic functions. 
The hydrodynamical interpretation of Hamel’s theorem is as follows. For an incom- 


pressible fluid moving in the plane with vorticity w, we have 


u, +v, = 0, v, — U, = 2w, (2) 


and there is a stream-function y such that 
u=—y,, v=y., Ay = 2w. (3) 


Thus Hamel’s theorem tells us that in order that a given distribution of vorticity may be 
consistent with vanishing velocity on the boundary B (the usual boundary condition 
for a viscous fluid in a fixed container), it is necessary and sufficient that 


[ov dx dy = 0, (4) 


U being an arbitrary harmonic function. 

However, inspection of Hamel’s proof (loc. cit. p. 266) shows that he made use of a 
Green’s function of the second type, i.e. a harmonic function G, with a singularity log r 
and making dG,/dn = 0 on B. There is, of course, no such function for Laplace’s equation, 
since this singularity and this boundary condition are inconsistent. 

Not knowing of Hamel’s work, I obtained Hamel’s result in 1935 in a rather special 
case (Proc. London Math. Soc. 40 (1935), 23-36) in a different way.** In the present note 
the theorem is extended to include compressibility. 

Theorem: A compressible viscous fluid moves inside a fixed connected boundary B, on 
which the velocity vanishes. An expansion 6(z,y) and a vorticity w(x,y) are consistent 


with this boundary condition if, and only if, 


[ ov - Vided <6, (5) 


*Received November 16, 1949. 
**Footnote added in proof (Feb. 20, 1950): The result (4) has recently been proved by J. Kampé de 
Fériet (Math. Mag. 21, 71-79(1947); Ann. Soc. Sci. Bruxelles (I) 62, 11-18(1948). 
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where U is an arbitrary harmonic function and V the conjugate harmonic function, such 
that 
U,=V,, U,= —-V.. (6) 
In purely mathematical language, equation (5) is a necessary and sufficient condition 
for the consistency of the equations 
U, tv, = 8, Vv, — U, = 2a, (u)z = 0, (v)p = 0. (7) 
Proof: Let 1, m be the direction cosines of the outward normal to B. Let 6 and w be arbi- 
trarily assigned. Let wu’, v’ satisfy 
ur+ov; = 8, vi — uy = 2, (lu’ + mv’), = 0. (8) 
It is well known that the solution (u’, v’) is unique, since the two partial differential 
equations define (u’, v’) to within the gradient of a harmonic function, and the normal 
derivative of the latter on B is then given by the last of (8). 
Denoting the integral in (5) by J, we have 


Il 


/ (QwU — 6V) dx dy 


I 


/ [U(v, — ui) — V(ul + v})] dx dy 


(9) 
és / [u’(U, + V.) + o'(V, — U,)] dx dy 
+ [ [U(l’ — mu’) — V(lu’ + mv’)] ds 
“B 
or, by (6) and (8), 
[= | U(l’ — mu’) ds. (10) 
“B 
Now if J = 0 for arbitrary harmonic U, it follows that 

(lv’ — mu’), = 0,7 (t1) 


since the values of U on B may be arbitrarily assigned. Combining (8) with (11) we get 
(u’)» = 0, (v’)p = 0; thus the condition I = 0 is sufficient. 

On the other hand, if (7) are consistent, then (u’)p = 0, (v’)g = 0, and so, by (10), 
I = 0; thus the condition I = 0 ts necessary. 

We get Hamel’s theorem on putting @ = 0. 
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BOOK REVIEWS 


On the theory cf stochastic processes and their application to the theory of cosmic radiation. 
By Niels Arley. John Wiley & Sons, Inc., New York, 1948. 240 pp. $5.00. 


This book aims at extending the general theory of stochastic processes, as developed by Kolmogoroff, 
Feller, Lundberg and others, and applying the resulting generalization to the theory of cosmic ray 
showers. The work is thus naturally divided into two parts, which as the author says in the Preface, 
may be read independently of each other. 

The mathematical section contains an elegant and rigorous, albeit somewhat repetitious, discussion 
of stochastic processes,—processes involving one or more random or statistical variables, the probability 
distribution of which depends on one continuously varying parameter. Although reference is made for the 
sake of generality to events labeled by Borel point sets in Euclidean space with a finite number of dimen- 
sions, the theory is immediately specialized to so-called Markoff chains, in which the stochastic variables 
take on only an enumerable manifold of values. 

The central problem in the theory is the existence of solutions of the diffusion equations for the 
required probability distributions, given the elementary transition probabilities. Arley discusses this 
problem in a quite general way, using a convenient matrix representation, in which the row and column 
labeling represents respectively the final and initial states of the stochastic variables. He introduces the 
concept of absolute exponentiability, namely the existence of a matrix 


© 
>> K"(t — s)"/m! 
m=0 
where K is a bounding matrix of the fundamental transition probability operator. It is then shown, except 
for normalization, that solutions of the diffusion equations exist for absolutely exponentiable types of 
operators. One of the contributions of Arley is to introduce a new set of sufficient conditions for absolute 
exponentiability, in addition to conditions established by Feller, Kolmogoroff, and Lundberg, whose 
results are included in detail in this monograph. Arley’s conditions are, indeed, just suitable for the cosmic 
ray problem. 

The problem of proving the normalization of the fundamental solution, (which was obtained using 
Peano-type product-integrals) is a difficult one, and Arley succeeds only in proving it for his “generalized 
elementary processes’’, namely, processes for which in each upward change, the stochastic variable cannot 
change by more than a certain number / of units. 

This section of the book will be valuable to statisticians and mathematicians interested in extending 
the general existence theory for stochastic processes. It is not of much practical use for applications, 
since as the author himself says, the existence of solutions is generally assured in physical problems, and 
the equations may be written down from physical considerations. Although, as stated on p. 46, the 
product-integral is the ideal mathematical tool for the theory of discontinuous processes, it is rarely of 
help in numerical work. 

The author introduces the second part of the work by saying (p. 15): “Jt is, namely, a fundamental, 
empirical law of nature that the region of validity of the simplifying approximations of which any theory must 
necessarily make use, is always far wider than might be justified by theoretical arguments.”’ (Author’s italics). 
The simplifying approximations made here are indeed extensive, and while the agreement with experi- 
ment is good, and perhaps more than would be expected at first glance, careful discussion would reveal 
sufficiently valid reasons for the agreement to make such a bold assertion unnecessary. In any case, it 
could stand further discussion. 

The second part of the work deals with cosmic ray showers, which are produced by the combination 
of the radiation of high-energy photons by electrons, and the production of electron pairs by photons, in 
their passage through matter, with the atomic collisions by the electrons accounting for the loss of energy 
and ultimate disappearance of the shower. A brief but adequate description of this process is given in 
the text. The most serious simplification made in the theory is the disregard of the energy-dependence 
of the various probabilities involved. The inclusion of this dependence would throw the problem outside 
the realm of application of Arley’s theory, and require a further generalization. 
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The author is aware of this difficulty, and proposes to overcome it by making a model in which 


particles duplicate themselves with a constant probability, but ‘die’? with a probability proportional 


to their depth in the shower-evoking material. It is hard to see the extent of the error produced by using 

this model instead of one in which particles ‘‘die’”’ when their energy is reduced to zero. 
Another simplifying assumption is to make the two duplication processes symmetrical 

in order to make the equations 


electrons 
each producing two photons only, to correspond with the reverse process 
amenable. Even then it is only possible to find the first and second moments of the desired distribution, 
which is the probability of observing n electrons at depth t. So a “reasonable” distribution with two adjust- 


able parameters is used, the Polya distribution 


b "/—1/b 
sin) = n y( | ‘Ja 4 1/b 
e\l + bn nv 


and for each value of t, b and # are chosen to yield the calculated first and second moments, fi and #2, 
The reasonableness of this procedure is shown, at least for large i and small t, by comparison with a 

solution obtained by numerical integration of the basic diffusion equations. The theory is applied to the 

by using the fairly well-established calculations of Bhabha and Heitler for the 


actual cosmic-ray case, 
death-rate’’ of the model 


latter case, of the first-moments problem (the average numbers), adjusting the 
to yield the same value of the number of particles at the maximum, then using the model to obtain the 


“c 


second moments, and finally using the Polya distribution as indicated. 

The results are first compared with Poisson distributions, which clearly are not valid owing to the 
correlated nature of the statistics, and are applied to the Geiger-counter shower measurements according 
to Rossi. The probabilities have to be integrated over an assumed initial spectrum of photons and elec- 
trons, which further complicates the calculation and makes difficult an evaluation of the validity of the 
procedure. The comparison with experiment is, however, fairly good. 

It should be stated, parenthetically, that much better average-number calculations for shower theory 
are now available (H. 8. Snyder, Phys. Rev. 76, 1563 (1949)), and that in particular the low energy calcu- 
lations of Arley and Ericksen may be off by a large factor. 

The general result is the demonstration of a practical way of predicting, albeit roughly, the shape of 
the Rossi curves and the results of similar statistical experiments. Little light, however, seems to have 
been thrown on the general and difficult analytical problem of the fluctuations in cosmic ray showers. 
It has been shown (Scott and Uhlenbeck, Phys. Rev. 62, 497, (1942)) that the second moments and the 
fluctuation can be calculated directly in the actual case, and an approach of this sort may prove more 
helpful in the general theory. 

The typography of the book is excellent, and there seem to be no essential misprints. The English 
is good, giving one but few reminders that the author is Danish. The photo-offset reproduction has 


reduced the original about 10% with little loss of clarity. 
C 
W. T. Scorr 


Mathematics dictionary. By A. A. Alchian, E. F. Beckenbach, C. Bell, H. V. Craig, 
G. James, R. C. James, A. D. Michal, I. 8. Sokolnikoff. Edited by Glenn James and 
Robert C. James. D. Van Nostrand Company, Inc., New York, 1949. v + 432 pp. 
$7.50. 


This work is based on the ‘James Mathematics Dictionary” which was devoted to the terms of 
elementary mathematics. To those have now been added the basic terms in “metric differential geo- 
metry, theory of functions of real and complex variables, advanced calculus, differential equations, 
theory of groups, theory of metrices, theory of summability, point-set topology, general analysis (re- 
ferring in particular to integral equations and the calculus of variations, analytic mechanics, theory 


of potential, and statistics.” 
W. PRAGER 
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Description of a relay calculator. By The Staff of the Computation Laboratory. Harvard 
University Press, Cambridge, Massachusetts, 1949. xvi + 366 pp. $8.00. 


This book gives a complete description of the well-known Mark II calculator which was built for 
installation at the Naval Proving Ground, Dahlgren, Virginia by the Harvard Computation Laboratory 
under the direction of Professor H. N. Aiken. 

The machine, now in operation at Dahlgren, is a general purpose digital computer with an inner 
relay memory of 96 registers together with 24 sets of switches for the storage of constants, and 2 function 
tables for evaluating certain elementary functions, together with 4 interpolator mechanisms; it has 2 
adding units and 4 multiplying units and it receives its orders from a punched tape. It has been so 
designed that it is possible to split it into two machines and operate each independently of the other. 

The code and control are so arranged that in each second the calculator can form 2 multiplications, 
1 additions and 6 numeric transfers. 

The body of the book is devoted to a careful and detailed description of these various units and of 
their modes of operation. It is clear that anyone who is interested in carrying out a calculation on the 
machine will find this work indispensable. It should also prove interesting to all those engaged in machine 
design since some of the points of view expressed are of importance in all branches of the design problem. 

The last chapter of the book is devoted to a discussion of problem preparation. After the statement 
of some general principles a number of examples are considered in detail; these serve to give the reader a 
grasp of the system of coding used. 

Herman H. GoupsTINE 


Extrapolation, interpolation, and smoothing of stationary time series with engineering 
applications. By Norbert Wiener. The Technology Press of The Massachsuetts 
Institute of Technology, John Wiley & Sons, Inc., New York and Chapman & Hall, 
Ltd., London. ix + 163 pp. $4.00. 


This monograph is the second of a series. The first, Cybernetics (Wiley, N.Y., 1948; Math. Rev. 9, 
598), presented the philosophy of “control and communication in the animal and the machine’’. This book 
presents in detail some of the methods and techniques for cybernetics. In particular, the book considers 
problems of statistical prediction and filtering with applications to the design of communication systems. 
It was first published as a classified NDRC report and has been supplemented with two appendices by 
N. Levinson. The author was the first to formulate problems in communication theory on a statistical 
basis and his methods are the alloying of methods from the two fields: communication engineering and the 
statistical study of time series. 

After introducing general concepts and summarizing some basic mathematical notions, the 
author continues in a readable and orderly fashion to explain the method. He introduces numerous exam- 
ples to illustrate its application. The first problem studied is that of linear prediction for single stationary 
time series. The past of a message up to and including time ¢ and certain statistical information about the 
message are assumed known; the message is a time series and the statistical parameter is its autocorrela- 
tion coefficient. He then shows how to select from a class of physically realizabie operations on the mes- 
sage that operation which gives the best prediction of the message at time ¢ + a@ in the future. The 
allowable operators are linear, invariant with respect to the origin of time, and the result of the operation 
depends only upon the past of the message. The best prediction is in the sense of least mean (over all t) 
square error. The best operator depends only upon the autocorrelation coefficient; it is best for the whole 
class of messages which are statistically alike in that they have the same autocorrelation coefficient. 
Noise is then assumed to be superimposed on the message. In a similar manner a method is given for 
finding the linear operator which gives the best approximation of the message at time ¢ + a. Here the 
best operator depends only upon the auto-correlation and cross-correlation coefficients of the noise and 
the message. The degree of lag (a negative) which is advantageous in filtering is discussed; a positive 
corresponds to the combined problem of filtering and predicting. The methods have the virtue that they 
do not ignore phase distortion and are not limited to steady-state behavior. The extension of the method 
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to multiple stationary time series is indicated. Application of the method to approximate differentiation, 
interpolation and extrapolation are discussed in the last chapter. Appendix A contains a table of Laguerre 
functions. Appendices A and B by N. Levinson discuss a simpler computational procedure and give a 
heuristic account of the author’s theory. 

Except for certain practical difficulties in computation and in the construction of circuits, the most 
obvious obstacle to the use of the theory is in making a choice of correlation coefficients suitable to the 
problem at hand. In this regard the engineer will need to gain more insight into the significance of these 
statistical quantities. Already this book has influenced the work of many people. Its public publication 
should widen the scope of its influence and should hasten the advancements which the author foresees 


as possibilities. 
J. P. LASALLE 

















